1 1 ,< *>s NfiVflL RESEARCH LOGISTICS QUflRTERLy JUNE 1970 VOL. 17, NO. 2 fe re 13 OCT 7 0% OFFICE OF NAVAL RESEARCH NAVSO P-1278 NAVAL RESEARCH LOGISTICS QUARTERLY EDITORS Rear Admiral H. E. Eccles, USN (Retired) The George Washington University F. D. Rigby Texas Technological College S. M. Selig Managing Editor Office of Naval Research Arlington, Va. 22217 O. Morgenstern Princeton University D. M. Gilford U.S. Office of Education ASSOCIATE EDITORS R. Bellman, RAND Corporation J. C. Busby, Jr., Captain, SC, USN (Retired) W. W. Cooper, Carnegie Mellon University J. G. Dean, Captain, SC, USN G. Dyer, Vice Admiral, USN (Retired) P. L. Folsom, Captain, USN (Retired) M A. Geisler, RNAD Corporation A. J. Hoffman, International Business Machines Corporation H. P. Jones, Commander, SC, USN (Retired) S. Karlin, Stanford University H. W. Kuhn, Princeton University J. Laderman, Office of Naval Research R. J. Lundegard, Office of Naval Research W. H. Marlow, The George Washington University B.J. McDonald, Office of Naval Research R. E. McShane, Vice Admiral, USN (Retired) W. F. Millson, Captain, SC, USN H. D. Moore, Captain, SC, USN (Retired) M. i. Rosenberg, Captain, USN (Retired) D. Rosenblatt, National Bureau of Standards J. V. Rosapepe, Commander, SC, USN (Retired) T. L. Saaty, U.S. Arms Control and Disarmament Agency E. K. Scofield, Captain, SC, USN (Retired) M. W. Shelly, University of Kansas J. R. Simpson, Office of Naval Research J. S. Skoczylas, Colonel, USMC S. R. Smith, Naval Research Laboratory H. Solomon, The George Washington University I. Stakgold, Northwestern University E. D. Stanley, Jr., Rear Admiral, USN (Retired) C. Stein, Jr., Captain, SC, USN (Retired) R. M. Thrall, University of Michigan T. C. Varley, Office of Naval Research C. B. Tompkins, University of California J. F. Tynan, Commander, SC, USN (Retired) J. D. Wilkes, Department of Defense, OASD (ISA) The Naval Research Logistics Quarterly is devoted to the dissemination of scientific information in logistics and will publish research and expository papers, including those in certain areas of mathematics, statistics, and economics, relevant to the over-all effort to improve the efficiency and effectiveness of logistics operations. Information for Contributors is indicated on inside back cover. The Naval Research Logistics Quarterly is published by the Office of Naval Research in the months of March, June, September, and December and can be purchased from the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402. Subscription Price: S550 a year in the U.S. and Canada, S7.00 elsewhere. Cost of individual issues may be obtained from the Superintendent of Documents. The views and opinions expressed in this quarterly are those ofthe authors and not necessarily those of the Office of Naval Research. Issuance of this periodical approved in accordance with Department of the Navy Publications and Printing Regulations, NAVEXOS P-35 Permission has been granted to use the copyrighted material appearing in this publication. MARKOV CHAIN ANALYSIS OF A SITUATION WHERE CANNIBALIZA- TION IS THE ONLY REPAIR ACTIVITY Alan J. Rolfe The Rand Corporation Santa Monica, California ABSTRACT This paper considers a group of S identical aircraft, each of which is partitioned into K parts which fail exponentially. The only way in which a failed aircraft can be repaired is by cannibalizing its out-of-commission parts from other failed aircraft. The evolution of the number of good aircraft over time is governed by the transient behavior of an absorbing Markov chain. We can therefore study this behavior by matrix multiplication although the computational problem grows large for K 5= 3. Some numerical results and some approxi- mations are also provided. 1. INTRODUCTION This paper presents a model of a system subject to cannibalization. The system is a simple special case of the very general coherent structure considered by Hirsch, Meisner, and Boll (HMB) in [3]. It differs, however, from the model of HMB in its treatment of time. In Sec. 2 the model is described in terms of a specific example. The reader can no doubt imagine a number of other situations where the same model would be appropriate. We also make clear how our model treats time differently from the model of HMB. Section 3 presents the Markov chain transition matrix which characterizes our problem, while Sec. 4 explores the computational chore associated with using the model. In Sec. 5 and 6, we exhibit several modest numerical examples and suggest means of lessening the computational chore for larger problems. Section 7 contains a special technique which is feasible for small-scale problems. Although this paper does not completely and explicity solve the problem it poses, we hope that the results given will provide quantitative insights into the effects of cannibalization which will lead to further work in the field. 2. THE PROBLEM We consider a group of S identical aircraft operating under the following conditions: (a) Every aircraft includes K parts or assemblies. The time-to-failure of the ith part is exponentially distributed with distribution function, 1 — e~V. Futhermore, the K parts are assumed to be noncritical items so that no matter how many of the K assemblies fail during a mission it is always possible for an aircraft to return safely to base. (b) A failed part can neither be repaired nor resupplied. However, there is sufficient repair capa- bility to cannibalize unfailed parts from other failed aircraft in an attempt to minimize the number of out-of-commission aircraft. We assume that initial supplies of spare parts are zero, although this condition can be relaxed as we indicate later. (c) The aircraft fly missions of fixed length T in which every "good" aircraft flies every mission. We assume that there is enough time between consecutive mission takeoff times for all aircraft to 151 152 A. J. ROLFE return from the previous mission and for all cannibalizations to be performed. We indicate below how the "every good aircraft flies every mission" rule can be modified. The basic quantity which the analysis below will provide is the expected number of good aircraft after t missions, t=\, 2, . . ., denoted by E[A{t)]. This information will, for example, allow us to measure the changes in force capability over time which cannibalization affords under the conditions above. It is important to distinguish the way in which our model differs in its treatment of time from that of HMB. In our case, only those parts which happen to be in the air during a given period of time are subject to risk. Good parts lying in partially cannibalized aircraft on the ground are not subject to failure. In effect, each part carries its own clock which is turned on when the part is sent up into the air and turned off when the part is on the ground. In the paper of HMB, however, every part that has not previously failed is subject to risk. There is a single clock which is turned on at time 0, and while it runs all parts are subject to failure. In Sec. 6 below we use the above fact to identify one of the results of HMB as a lower bound for E[A(t)]. 3. THE MARKOV TRANSITION MATRIX If we let rii(t) be the number of copies of part i which survive t missions, then n(t) = (ni(t), . . ., nA-(O) represents the number of parts of each type which survive t missions. This is the basic state vector of our process. As n.i{t) is an integer and ^ rii{t) ^ S for i=l, 2, . . ., K, the process has (S + 1) K states for each t. However, for our purposes we can lump all states having one or more ni{t) =0 into a single "zero" state, thereby reducing the number of states to S K + 1. It is readily seen that the process is Markovian and that its evolution from mission t to t+l can be described by a one-step Markov chain transition matrix P of dimension S K + 1 by S K + 1. If the rows and columns are ordered lexicographically, P will be a lower triangular matrix with 0's above the principal diagonal. It is convenient to label row and column m by d m , where d\ =(zero state); <k =(1, . . ., 1); d S K +1 =(S, . . .,S). Then the matrix element py, where dt= (t'i, . . ., i/<), dj= (j\ . . .,Jk), and M,= min ii , is given by i=\ K / M. \ (2a) f[ ' ,„ ) (e- K i T ) j i- i ' +M i(l-e- K i T ) i i-Ji i jJi \ji — ii + Mj/ if i,y>l and;, =S i, <;, + M„ for all /=1, . . .,K; (2b) if i, j > 1 but ii < ji or i t > ji + Mi for at least one / = 1 , . . . , K; (2c) lifi=j = l; (2d) Oif £=1,j>1; MARKOV CHAIN CANNIBALIZATION ANALYSIS 153 i (2e) 1~2 Pidi£i>l,j=l- rf=2 These entries become apparent if we realize that Mi copies of each part are tested on each mission and that each part is assumed to fail independently of other parts. Once the zero state is reached the "experiment" ends as all aircraft are out-of-commission. 4. THE COMPUTATIONAL CHORE If we let ri(<) be the row vector of unconditional state probabilities after t missions, then from Markov chain theory we know that (3) n(*)=n(0)F, «=i,2, . . ., where 11(0) = (0, . . .,0,1) as the process starts in state (S, S. . . . , S) and the states of II() are ordered lexicographically as are those of P. Furthermore, P represents an absorbing Markov chain [4] with d\ the absorbing state and all other states transient. From (3) we see that to go from Tl(t) to Yl(t-\-\) we must multiply 11(0 , a 1 byS K + 1 matrix, byP, anS^-f 1 byS A + 1 matrix. This operation re- quires roughly (S 2K /2) multiplications and the same number of additions. The major problem in using the model is one of computational feasibility as (S 2A /2) can be a large number. For example, if S = 24 and K = 3, (S 2K /2)= 1.5 X 10 8 . Thus it appears that either S or K must be fairly severely restricted if we are to perform the computation in an acceptable amount of time; the following sections describe some numerical results which suggest ways of lessening the computational chore. The first two moments of the number of good aircraft after t missions are given by (4) E[A(t)]= S 2 nj(t)Mj and j=1 (5) E[A(t)*] = S! 2 Uj(t)Ml where Mj= min {ji\dj={ji, . . ., >*)}• 1=1, . . .,K Since Mj is independent of t for each;, the computation of E[A{t)] should proceed rapidly once Yl(t) has been determined. As the sole absorbing state d\ of the Markov process has Mi = it is clear from (4) that lim E[A(t)]=0; this is of course what one would expect given the no-repair situation being modeled. The computation of 11(f) would generally be iterated until E[A(t) ] was fairly close to zero. 5. NUMERICAL RESULTS FOR K = 2 In order to simplify the computational task, we first limit ourselves to the case of two cannibalizable parts per aircraft. There remain, then, two dimensions to the problem: S, the number of aircraft under test, and g = {e~ x,T , e~ k2T ), the vector whose components are the probabilities of nonfailure of each part on one mission of length T. Table 1 presents results for q fixed and S variable, while in Table 2 q is variable and S fixed. In both cases the quantity of interest is E[A(t)], the expected number of good aircraft after t missions. 154 A. J. ROLFE Table 1 \ s t \. Any S, no cannibalization; (0.72)' 2 5 10 13 1 0.720 0.734 0.762 0.777 0.783 2 0.518 0.557 0.600 0.622 0.628 3 0.373 0.431 0.480 0.500 0.504 4 0.268 0.336 0.384 0.402 0.405 5 0.193 0.263 0.308 0.323 0.326 6 0.139 0.206 0.247 0.260 0.261 7 0.100 0.161 0.198 0.209 0.210 8 0.072 0.125 0.159 0.168 0.169 9 0.054 0.097 0.128 0.135 0.136 10 0.037 0.075 0.102 0.108 0.109 11 0.027 0.058 0.082 0.087 0.087 12 0.019 0.044 0.066 0.070 0.070 13 0.014 0.034 0.053 0.056 0.056 14 0.010 0.026 0.042 0.045 0.045 15 0.007 0.019 0.033 0.036 0.036 (Entries are E[A(t)]IS; q = (0.9, 0.8). For the no cannibalization case E[A(t)]/S = (0.72)' results from the fact that 0.72 is the probability that one aircraft survives one mission and that we have a sequence of t bi- nomial trials.) From Table 1 we see that for a fixed value of t, the ability to cannibalize produces an increase in E[A(t)]IS which also increases with S. However, we note that as S increases the quantity E[A(t)]/S Table 2 t \v 0.9, 0.8 0.8, 0.8 0.7, 0.7 0.6, 0.6 0.8, 0.5 0.5, 0.5 0.9, 0.4 0.5, 0.4 1 2 3 4 5 6 8 10 12 14 0.777 0.804 0.803 0.806 0.805 0.802 0.803 0.804 0.804 0.804 0.731 0.768 0.773 0.775 0.777 0.776 0.776 0.775 0.773 0.620 0.661 0.667 0.668 0.668 0.667 0.660 0.514 0.558 0.563 0.558 0.556 0.550 0.494 0.506 0.504 0.503 0.503 0.412 0.456 0.456 0.450 0.399 0.400 0.401 0.400 0.353 0.391 0.389 0.382 (Entries are E[A(t)]/E[A(t- 1)]; S = 10.) converges for values of S around 10. The entries for S — 13 are not much different from those forS = 10. Hence it seems that we can learn as much about E[A(t)]/S by studying 13 aircraft as by studying 24 or 48 and the computational job is markedly easier. Of course, this conclusion is strictly valid only for K — 2 and q = (0.9, 0.8). A question of interest (as yet unanswered) is, "At which value of S does con- vergence occur, given a particular K and q?" In Table 2 we exhibit the behavior of successive values of E[A(t)] for a group of 10 aircraft and a variety of g-vectors. One fact is immediately apparent: For t 2= 2, E[A(t)]/E[A(t— 1)] =p, p a constant. MARKOV CHAIN CANNIBALIZATION ANALYSIS 155 This is equivalent to saying £'[/i(f)]/S = p'" 1 , t Ss 2, which implies (6) E[A(t)]IS = p ( P y-\t^l, where pomEtA(l)]IS. This provides an immediate analogy to the no-cannibalization case, where we know that E[A{t)]/S = A.', «Sl, with A = <7i<7 2 where q = (qi , q 2 ) . (It was found that for S = 5 the geometric behavior was less stable than that characterized by Table 2 for S= 10). The next section will suggest several possible ways around the computational problem. 6. APPROXIMATIONS FOR K>2 In Sec. 5 we indicated that for K= 2, any S 3= 10 and any q we could approximate the behavior of E[A(t)]/S by studying that of the case S— 10. A typical computer run on the 360/65 for S= 10 and q= (0.9, 0.8) took approximately 7 sec. If K were increased to 3, the information in Sec. 4 would suggest a run time of roughly 10 min., still a feasible amount of time. Hence for K= 2,3 we can study a group of 10 aircraft and be fairly confident that the E[A (t)]IS so obtained closely approximates that for any larger value of S and the same q. For K 25 4, however, studying even 10 aircraft becomes prohibitively costly as each time K is increased by 1, computation times are increased by a factor of 100. Thus, although the geometric behavior exhibited by E[A(t)]/S in (6) would probably still hold, the problem of determining p and p could no longer be handled by numerical means. An approximation for E[A(t)] is given by (7), with q= {q u ■ ■ • , <7a ) and a=mjn </,, (7, 5f*^,„i Relation (7) is based on the imprecise reasoning that A(t) is equal to the number of survivors of the part-type for which qi is minimized. The numbers in Table 2 suggest that the RHS of (7) is often an upper bound for E[A(t)]IS rather than an approximation; however, column 5 of that table shows that this is not always so. For t=l, the relationship E[A( 1)]/S =S a follows readily from min Xj^Xi, i=l, . . . , K, i = 1 K K A lower bound for E[A(t)]/S is given by (8), with |8 = r[ qi, The referee suggested the following better lower bound, which can be derived from formulae 4.3 and 4.9 of [3]. This result can also be derived directly without reference to [3] and is based on the supposition that every surviving part be tested on every mission. (9) E[A(t)] The referee provided the results shown in Table 3 whose entries are the lower bound, (9), for the case shown in Table 1. Expression (9) appears to be an excellent lower bound in this instance. 156 A. J. ROLFE Table 3 \ s t \^ 2 5 10 13 1 0.734 0.759 0.777 0.782 2 0.554 0.598 0.620 0.626 3 0.423 0.473 ■ 0.495 0.501 4 0.323 0.374 0.396 0.401 5 0.247 0.296 0.316 0.320 6 0.187 0.233 0.252 0.256 7 0.142 0.184 0.201 0.204 8 0.106 0.144 0.160 0.163 9 0.080 0.113 0.127 0.130 10 0.059 0.088 0.101 0.103 11 0.044 0.068 0.080 0.082 12 0.032 0.052 0.063 0.065 13 0.024 0.040 0.049 0.052 14 0.017 0.030 0.039 0.041 15 0.013 0.023 0.030 0.032 7. A RECURSIVE RELATION FOR THE GENERATING FUNCTION OF THE STATE PRORABILITIES The fact that the P-matrix, as given by expression (2), is lower triangular allows us to develop a recursive relation which, at least in theory, produces a closed form expression for the generating func- tion of the state probability vector, Yl{t). Inverting these generating functions presents a barrier if S and Kare large; however, we present the method as it might sometimes be desirable to study a problem with small S and K. Setting N=S K +l, we let {ptj, i, j=l, . . ., N} be the elements of the P-matrix. If II(t) = (IIi(f), . . ., Il/v(t)) is the vector of state probabilities, we define the generating functions (10) Ui(Z) = £ Yli(t)Z', i = l, 2, . . ., N. t = The basic Markov chain relation is (11) nu) = n{t-\)P. Using (10), (11), and the fact that pjj = 0,j> i, we can derive U N (Z) =p NN ZU N (Z) + 1 This is equivalent to (12a) (12b) Yl i (Z) = ^ Pji Zn j (Z), i=N-l,N-2, . . ., 1. n N {Z) = \j{\-p NN z) n,(Z) = (z 2 p ji n j (Z))i(i-p ji z), i=N-i,N-2, . . ., 1. j = i + l MARKOV CHAIN CANNIBALIZATION ANALYSIS 157 From (12a) we know at once that Us(t) = (p NN y, t^O. Then (12b) can be solved recursively as IT(Z) depends only upon flj(Z), j>i. The difficulty in a large problem is that inverting Yli(Z) is not trivial as it was in (12a) because the generating function will have a large number of multiple roots. There do exist both exact and approximate methods for such problems, some of which are found in [2, Chap. XI]. We present a short numerical example to illustrate the technique: S =2 K = 2 q= (0.5, 0.4) From expression (2), we find that P = 1 0.8 0.2 0.5 0.3 0.2 0.6 0.2 0.2 _0.52 0.24 0.08 0.12 0.04 Then expression (12) produces n 5 (Z) = i/u-o.04Z) fl 4 (Z) = (0.12Z)/(1 -0.2Z) (1 -0.04Z) n 3 (Z) = (0.8Z)/(1 -0.2Z) (1 -0.04Z) n 2 (Z) = (0.24Z)/(l-0.2Z) 2 (l-0.04Z) n I (Z) = (0.0016ZM5 + Z)(65-Z)/(l-Z)(l-0.2Z) 2 U-0.04Z). If we use methods in [1] to invert these generating functions and then use the relation we find E[AU)] = o{n l u)} + i{n 2 (t) + iu(t) + n 4 (t)} + 2{Ti i (t)}. ^(*)]=|«(.2V+|(.2)'+|(.04)',«s»0. Thus we have found a closed form expression for E[A(t)]. 8. VARIATIONS OF THE MODEL Several of the assumptions made in Sec. 2 can be relaxed without materially affecting the develop- ment outlined above. The only change is in the Markov transition matrix P: 158 A. J. ROLFE (a) There may be an initial supply, s,, of spares for part i, i= 1, 2, . . . , K. In this case the possible states of part i are 0, 1, . . .,5, + S and the P-matrix has f[( Si + S) + l rows and columns. This would make the computational task mentioned in Sec. 4 more difficult; (b) The assumption that "every good aircraft flies" every mission" can be replaced by the re- quirement that only Q aircraft, if available, fly every mission. This would be accomplished by sub- stituting Q for Mi in expression (2a) whenever Q < Mi. If this were done, the P-matrix would be made more sparse while E[A(t)] would approach zero less rapidly because the Markov process would take more time to work its way from the initial state to the absorbing state. 9. ACKNOWLEDGMENTS The author would like to thank Dr. D. M. Landi of The Rand Corporation for his part in bringing the model to the author's attention. The author is especially indebted to an anonymous referee who supplied many useful comments, all of which are incorporated in this paper. REFERENCES [1] D'Azzo,J. J. and C. H. Houpis, Feedback Control System Analysis and Synthesis (McGraw-Hill, New York, 1960). [2] Feller, W., An Introduction to Probability Theory and Its Applications (Wiley, New York, 1957), Vol. 1. [3] Hirsch, W., M. Meisner, and C. Boll, "Cannibalization in Multicomponent Systems and the Theory of Reliability," Na'v. Res. Log. Quart. 7.5, 335-359 (1968). [4] Kemeny, J. G. and J. L. Snell, Finite Markov Chains (Van Nostrand, Princeton, New Jersey, 1960). SOME ESTIMATES OF RELIABILITY USING INTERFERENCE THEORY M. Mazumdar Westinghouse Ktscarch Laboratories . Pittsburgh, Pennsylvania 15235 ABSTRACT According to "interference theory" of reliability, a component fails if the maximum stress exceeds the component's strength. Assuming that both these quantities are random and their distributions are normal, we obtain in this paper some point and interval estimates of reliability when the stress distribution is known and a few observations exist on com- ponent strengths. 1. INTRODUCTION In recent years, many writers have attempted to analyze the reliability of a component by the application of probability arguments to a certain physical model of failure. According to this model (which has come to be recently known in the literature as "interference theory"), a component fails if at any moment the applied stress (or load) is greater than its strength or resistance. The stress is a function of the environment to which the component is subjected, and its value at any point of time is considered to be a random variable. The strength of the component is measured by the stress re- quired to cause failure. Strength depends on material properties, manufacturing procedures, and so on. If the components under question are mass produced and their selection in a given system is assumed to be made at random, then the strength should also be considefed a random variable. The reliability of a component during a given period [0, T\ is taken to be the probability that its strength exceeds the maximum stress during the entire interval. The above formulation of the reliability problem was first proposed by Birnbaum [1], who also gave a method for estimating the related probability. In recent years, the model has found an increasing number of applications in many different areas, especially in the structural and aircraft industries. The stress to which a component is subject in its operating environment can often be predicted using the available technological knowledge about the mechanical, thermal and other relevant condi- tions of the system, and the way they interact with each other. In many cases mathematical models are available which allow one to compute deterministically the values of "stress" at different points of time corresponding to a given set of input values describing the initial conditions of the system. For example, as Church and Harris [2] have pointed out, in a missile flight these initial values may cor- respond to the propulsive force, angles of elevation, changes in atmospheric conditions, etc. Since these input values have statistical distributions associated with them, the calculated "stress" values will also be randomly distributed. Moreover, since it will be possible to generate a practically unlimited number of such stress values using simulation methods by varying the input values in accordance with their respective distributions and hence to reconstruct the stress distribution as accurately as desired, we might assume that the distribution of the maximum stress is known. The strength distribution can- not, however, usually be computed from such a priori considerations. It can only be estimated using 159 160 M. MAZUMDAR statistical methods from the results of tests carried out for this purpose. The number of such tests is likely to be small in most cases due to practical limitations. Under, the assumption that the distributions for strength and the maximum stress are both normal, we derive in this paper some point and interval estimates of reliability using "interference theory," when as above the stress distribution is known and the strength distribution is unknown, but a few observations on the latter are available. The point estimates are given in section 2, and the interval estimates are derived in section 3. Owen, Craswell, and Hanson [7] considered a similar problem and gave confidence intervals for the reliability when both stress and strength are (a) jointly bivariate normally distributed and observa- tions are in pairs, or (b) when they are independently normally distributed, but they have a common (but unknown) variance. Most recently, Church and Harris [2] considered the same problem and provided a procedure for obtaining confidence intervals, which they compare with that of Govindarajalu [3]. Their procedure, although empirically demonstrated to be superior to that of [3], is, however, inexact since it uses the asymptotic normal approximation of a given statistic and requires the substi- tution of the population mean and standard deviations by their observed sample values. Our procedure for obtaining the confidence interval is based on that given in [7] and, though exact in a certain fashion, is admittedly inefficient. Some comparisons of this procedure with that of Church and Harris [2] are given in section 3. Sampling experiments show that the average length of the interval associated with the Church and Harris procedure is smaller than that of our procedure although it appears that, in general, it fails to cover the true parameter value more often. 2. POINT ESTIMATES Let X denote the maximum stress to which the component is subject over a given period of time. We assume that X is normally distributed with mean /jli and variance cr'f, and that both these values are known. Let Y denote its strength. We assume that Y is also normally distributed with unknown mean /x-z and variance erf. In what follows, we shall first assume that cr| is known and then that it is unknown. No doubt the latter is the more realistic case. In terms of our definitions we readily see the reliability to be (1) p = Pr{Y-X>0} = <b( ^-^ ) \ Vo-f + o-f/ ' where ®(t)=—^= (' e- x *' 2 dx. Let Yi,Yz, . . . , Y m denote m mutually independent observations on the strength of the components. These values are obtained from bench tests. On the basis of these observations, we obtain below the minimum variance unbiased (MVU) and the maximum likelihood estimates of p. 2.1 Case 1: cr| known Let Z be a normal random variable with mean and variance erf, and let (2) / = 1 if Y x + Z> /x, if Yi + Z^m RELIABILITY BY USING INTERFERENCE THEORY 161 Then (3) E{I}=Pr{Y l +Z> f x 1 } =P- It therefore follows from Blackwell-Rao and Lehmann-Scheffe theorems (See Chap. 3. [5]) that _ _ »' the conditional expectation E{I\Y} is the MVU-estimate of p, where F=VF;/m is the sufficient i= 1 statistic. Since moreover the conditional distribution of Y t + Z given F is normal with (4) £{F,+Z|F} = F, Var{F 1 +Z|F} = cr2 + o-|('l--) , the MVU-estimate of p in this case is given by (5) p, = E{h\Y) = Pr{Y i + Z>p. 1 \Y} = *( ^ ). VVo-f + o-fd-l/m)/ The maximum likelihood estimate of p, on the other hand, is given by It is seen that as m becomes large, the two estimates become equivalent. 2.2 Case 2: erf unknown In this case the sufficient statistics is the pair (F, S.J) where »' _ Si=2 (F,-Y)*. Therefore, using arguments identical to those of Case 1, it follows that in this case, E{I/Y, S'i} is the MVU-estimate of p. From [6] we find the conditional density of Y\ given F and Sf to be ' m \ 2 ) 1 f (Y-yV m ] 2 m — ,- „/m — 2\ S-» 1 \ S-> J m — otherwise. 162 M. MAZUMDAR Therefore the MVU-estimate of p in this case is (8) p n = E{I\Y,Si} = Pr{ Y x +Z >/xi | Y,S%} V m I i Y—y / m By use of the transformation 2 = - + - — - — \ - the integral on the right hand side of (8) reduces to z z 02 v m — I Jo r (9) />// = I ' . r( o m , 2) . zf"-'"/ 2 (l-z) (, "- 4)/2 ^(« + j82)ffe, 'm — Z\ v (m — Z x where a =(s.,J— 1- Y—pn )/cr,, m and /3 : 2S 2 \m — <j\ \ m It does not appear that this integral can be evaluated in a closed form; however, it can be evalu- ated very easily for particular examples by the use of numerical integration. Thf* maximum likelihood estimate forp in this case is (10) P// = <*>'- where Vaf+4 S? ss = 2 m-r pu is biased, but might be quite efficient when m is reasonably large. Also, it can be shown that for large values of m, the MVU-estimate pu and the maximum likelihood estimate/?// are equivalent. 3. INTERVAL ESTIMATES As in the discussions in the preceding sections, we give the interval estimates of reliability for two separate cases (a) when cr\ is known, and (b) when o" 2 is unknown. We derive here only the lower confidence limit because of its importance in the present application. An extension of our procedure to two-sided intervals will be obvious. 3.1 Case I: o~ 2 known Let K a be such that <$>(K a ) = a. Since Y is normal with mean p.-, and variance cr 2 /m, it follows that ,11) Pr{ V "<r-*> «-*.}-!-„, f fii — fiz Ui — Y—K a <T 2 lVm} . or rr \ , =£ \ = 1 — a, I V o-? + o- 2 V rr? + 0-2 J 'oi + 0- 2 RELIABILITY BY USING INTERFERENCE THEORY 163 Pr k <^1 h- Thus the 100 (1 — a) percent lower confidence limit for p in Case I is (12) p,,dl-a)=<t>l ./-—-; )• \ Vo-f + crj ' 3.2 Case II: &\ unknown The procedure suggested by Church and Harris [2] for obtaining the interval estimate in this case is based on the statistic v- ? -^ ^ai + sl They have shown that V is asymptotically normal with mean (ju, 2 — Pi )/ V^f^ and variance (/Lt 2 — (JLi) 2 <ri OV — m 2(m-l)(o-'f + o-2) 2 1 ■ "2> Using normal approximation for the distribution of V and upon replacing the unknown values of fx 2 and <t\ by their respective sample estimates Y ands*, they obtain (the symbol = denotes approximately equal to) (13) Pr { / *IZ^ > V+ K a CTv] = l-«, I Vo-2 + O"! J or Pr{p>4>(K+^„o-v-)} = 1-a, or Pr{p><D(r+^)}=l-a, where &l ai + ,\ ■+ m 2(m-\)((T 2 + siy Since this procedure involves two approximations both of which may not be close in the non-asymptotic case, the above results may not be satisfactory in small samples. We next suggest a procedure based on results given in section 4 of [7]. This procedure is "exact," but is admittedly inefficient because it does not make use of all the available information; however, the use of sampling adopted in this method seems to be the only way we can eliminate the nuisance param- eter &\ without introducing any approximations, and we believe that the properties of the procedure might therefore be of interest. It seems possible that ad hoc procedures can be developed which use this sampling concept, but are more efficient than the procedure described here. Results obtained from simulation experiments given below seem to indicate that our procedure covers, in general, the true parameter value more often than the Church and Harris procedure, although the average interval obtained from the latter procedure is significantly shorter. According to our procedure, we first draw, using random number tables, a sample of m independent observations form a normal distribution with mean and variance cr\. Having taken these m observations 164 M. MAZUMDAR we consider the transformed observations It is clear that (14) and Now let (15) and Then (16) Wi = Yi + Zi, i = l, 2, . . ., to. E{Wi}=p 2 Vai{Wi}=(r\ + trl w=^ — TO £ {Wi-wy m^\ Sw Sw Vof + o-f Vo-, 2 + o- 2 2 J/ VqHd/ is distributed as a non-central < with degrees of freedom = to— 1 and the non-centrality parameter 8— Vto(ju,i — /A2)/ Vo-f* + o"| . Let the left hand side of (16) be denoted by T,,, 1. It then follows from the formulas given in [4] that (17) Pr{ ^ (/"-**»> >S(m-l < r ffl - t , q) j = l-q, ( Vo- 2 + crl J where from Eq. (33) of [4], 8(to— 1, T m -u ot) is defined as (18) 1 8(to — 1, r m _i, a) =T m -i — \(m— 1, r m _i, a) 1 + 71- , 2(to-1) 1/2 and k(m— 1, r m _i, a) can be obtained from Table IV of Johnson and Welch [4]. From (17) and (18) it now follows that 5(to — 1, T m -u a) (19) Pr {p>^( TO = l-a, so that a lower (1 — a) percent limit forp following this procedure is given by <J>(S(to— 1, T m - U a)/Vm). RELIABILITY BY USING INTERFERENCE THEORY 165 We now exhibit the results obtained from 200 artificial sampling experiments comparing the above two procedures using the same set of random numbers. In Table 1 the Church and Harris procedure is denoted as Procedure A and the method suggested in this paper as Procedure B. The particular sample sizes shown in the table were chosen because of the comparative ease with which the values of 8(m — 1, T m -i, a) can be evaluated from the existing tables. Table 1. Comparison of Confidence Interval Procedures 0*i= 0, o-, = l; 1 -a = 0.95) p p-z 0"2 m Number of one-sided con- fidence intervals out of 200 which covered true param- eter value Average length of 200 confidence intervals Proc. A Proc. B Proc. A Proc. B 0.90 1.43287 0.5 10 184 191 0.15 0.28 37 183 195 0.12 0.19 1.81246 1.0 10 171 187 0.19 0.26 37 185 193 0.14 0.19 4.05278 3.0 10 187 183 0.25 0.27 37 187 191 0.17 0.18 0.95 1.83905 0.5 10 181 188 0.087 0.19 37 188 190 0.067 0.10 2.32624 1.0 10 184 194 0.12 0.21 37 184 190 0.08 0.11 5.20163 3.0 10 187 189 0.18 0.20 37 185 189 0.10 0.11 0.99 2.60055 0.5 10 178 187 0.022 0.095 37 186 193 0.015 0.038 3.28946 1.0 10 182 194 0.045 0.097 37 183 189 0.022 0.039 7.35546 3.0 10 190 193 0.093 0.11 37 192 188 0.034 0.039 REFERENCES [1] Birnbaum, Z. W., "On a Use of the Mann-Whitney Statistic," Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability (Uni\ersity of California Press, Berkeley, Calif., 1956), Vol. 1, pp. 13-17. [2] Church, J. D. and Bernard Harris, "The Estimation of Reliability from Stress-Strength Relationships," Technometrics, 12, 49-54 (1970). [3] Govindarajalu, Z., "Distribution-Free Confidence Bounds for P{X < Y}," Annals of the Institute of Statistical Mathematics, 20, 229-238 (1968). [4] Johnson, N. L. and B. L. Welch, "Applications of the Non-Central t Distribution," Biometrika, Vol. .31, 362-389 (1940). [5] Lehmann, E., Notes on the Theory of Estimation, Mimeo Notes (University of California Press, Berkeley, Calif.) (Reprinted 1962). [6] Lieberman, C. J. and G. J. Resnikoff, "Sampling Plans for Inspection by Variables," Journal of the American Statistical Association, 50, 457-516 (1955). [7] Owen, D. B., K. J. Craswell, and D. L. Hanson, "Non-Parametric Upper Confidence Bounds for P{Y < X) and Confidence Limits for P{Y <X} When X and Y Are Normal," Journal of the American Statistical Association, 59, 906-924 (1964). THE OPTIMAL BURN-IN TESTING OF REPAIRABLE EQUIPMENT* John M. Cozzolino Department of Statistics and Operations Research University of Pennsylvania Philadelphia, Pennsylvania ABSTRACT The subject of this paper is the utilization of the "infant mortality" or decreasing failure rate effect to improve the reliability of repairable devices. Decreasing failure rate implies the possibility that devices which exhibit it can be improved by "burn-in testing" of each unit. Such a test serves to accumulate operating time while shielded from the full costs and con- sequences of failure. A general formulation of the burn-in test decision for repairable devices is presented and some special cases are solved. A class of models, indexed by the degree of partial replacement present in the repair process, is considered and numerical results for the optimal policy are given for several members of that class. A comparison of those results reveals the profitability of testing increases with the complexity of the repairable device. 1. INTRODUCTION The subject of this paper is the utilization of the "infant mortality" effect to improve the reliability of repairable devices. The term infant mortality denotes a decrease with age of the conditional proba- bility of failure of a device and, consequently, an improvement in reliability with time as the device operates. This phenomenon is widely observed. For example, decreasing failure rate has been found to occur in electronic systems and components [12], biological systems [5], social systems such as business enterprises [20], and also, in mechanical systems. In the case of repairable equipments de- creasing failure rate is often observed after the occurrence of failure and repair. Decreasing failure rate implies the possibility that devices which exhibit it can be improved by "burn-in testing" of each unit. The term burn-in testing is taken to mean test stand operation under environmental conditions as close as possible to true operating conditions. The basic purpose of this test is to accumulate operating time while shielded from the full consequences and costs of failure. The measurement of reliability, while not a primary goal, may be a partial by-product. The main theme of this paper is a general formulation of the burn-in test decision problem for repairable equipment. A continuous stage dynamic program is used to treat this sequential decision problem and some solutions are given for special cases. Space does not permit a full discussion of the class of models for which this decision structure applies; however, one broad family of models, indexed by the degree of partial replacement present in the repair process, is presented. Numerical results for the optimal test policy are given for several members of that family, and a comparison of those results reveals an increase in the profitability of testing as a function of the complexity of the repairable device. *The work herein reported was supported by the Office of Naval Research under Contract Nonr 1841 (87); NR 042-230 while the author was at the Operations Research Center at the Massachusetts Institute of Technology. The computation reported here was done in part at the Computation Center at Massachusetts Institute of Technology. 167 168 J. M. COZZOLINO 2. LITERATURE REVIEW The distribution theoretic conditions under which the expected future life of a device increases with its age are shown by Watson and Wells [27]. Based upon the work of Knight [18], one stated result is that the expected future life of an aged device always exceeds that of a new one if the device has a monotonically decreasing failure rate. Several authors have investigated the burn-in test problem for the case of unrepairable equipments. In that problem a failed device is considered worthless. The quantities investigated are the burn-in time required to achieve a specified expected future life of the surviving units. Watson and Wells [27] discuss these quantities for several failure distributions. A nonparametric method is considered in a paper by Lawrence [19]. Assuming only a monotonically decreasing failure rate, he derives bounds upon the burn-in time required to achieve a specified expected future life. Recently, interest has developed in the burn-in testing of repairable devices. This potentially important problem requires some explicit assumptions about the nature of the relevant failure and repair processes. This necessity follows from the need to assess the relative worth of the repaired device. Kooharian and Young [28] consider a useful class of reliability models for repairable complex devices. The outstanding feature which they prove for this class is a special simplified form of the burn-in test decision policy. A similar model for the "debugging period" is given by Bazovsky [4]. Concern for the long run reliability improvement of complex systems is evidenced by a paper by Barlow and Scheuer [3] entitled "Reliability Growth During a Development and Testing Program." They attributed the improvement mainly to engineering changes. In a report upon which this paper is based [7], the author presents some models of decreasing failure rate processes incorporating explicit consideration of the repair process. This union of failure and repair models appears to be a fruitful basis for the discussion of burn-in testing in view of the interaction between the two processes. A complete picture of the life of the device is given in terms of a stochastic process, usually with non-independent interfailure times. 3. ELEMENTS OF THE PROBLEM The Unit Test A unit test starts with the testing of a new unit and ends with the output of one tested unit. If a failure occurs during test, the test is interrupted, the unit is repaired, and the test then resumes. The time elapsed during repair is not counted. The Failure Process for New Units Given a population of indistinguishable new devices, a specified and constant environment in which they are intended to operate, and a set of conditions which define failure, a distribution function Fo(tl6) with parameters 6 and density fo(t/6) is utilized t6 represent the uncertainty which exists re- garding t, the duration of operating time prior to the occurrence of failure for any one of the devices. The failure rate of that density, also known as the hazard rate, conditional probability of instantaneous failure, or force of mortality, is a very important function for determining the reliability characteristics of the device. For a distribution satisfying Fo(0~/6) —0, the failure rate, a function of time, is defined as It will be assumed throughout thatfo(t/6) is a continuous function of t and that Fo(O~/0) = 0. OPTIMAL BURN-IN TESTING « 169 Density of Future Time The original density, /o, represents the process only before testing begins. Suppose that at some time, r, after the test starts no failure has been observed in (0, t). Let t = future time to failure. Then the total time to failure is (£+ r). Of this total time to failure the only random part remaining is t. The relevant density function is that of future time to failure. This renormalized form of the density will be denoted here by f mr)- M,+Tle) - [l-F„(r/«)] While this density function usually represents a larger family of densities than does/oU/0), its failure rate is the same except displaced to the left by t; no(t + Tl8). The Repair Process We consider that a failure during test will be followed by a repair of the failed unit. The rule which assigns a new density of time to next failure for the repaired unit will be called the "repair rule." In many analysis it is assumed that the repaired unit has the same density of time to next failure as the new units. This practice ignores the information available from the failure and from the past history of the unit. Since it is possible to distinguish between new and repaired units we realize that it may be valuable to do so. This is particularly relevant for burn-in testing models since the purpose of testing is to create a difference between new and aged units. The only restriction which we shall place upon the new density is that it be from the same family of density functions as is f{t/8, t). Hence we allow t to have the meaning of a parameter of the density rather than strictly the time since last failure. The repair rule then is a transformation which determines which density of that family applies. The most general transformation to be considered here is one that allows the new parameter values to be a function of the parameter values which applied just prior to the failure. Hence the new parameter values are d R =o R (o,T) and tr = t h {6, t), where (8, r) are the parameter values just prior to failure. These functions will be required to have the property that, when applied to any point in the space of allowed parameter values, they will produce a point which is also inside that space. Also, their con- tinuity in t will be required. The repair rule results from an explicit model of the repair process. Examples of this are given in [71 The Reliability Value Function The reliability value function specifies the value that the tested unit would have if testing were terminated. It gives the value as a function of any possible test position (6, t). It will be denoted by: U(8, t) = the values of a unit when testing is terminated and a time when the state of process is {8, t). This class of functions, U(8, r), includes all expectations of the future behavior of the unit, since it includes all the parameters of the density function of future life. 170 J- M. COZZOLINO The Cost of Testing The last remaining element of the problem is the cost of testing. In this we seek to represent the actual cost of using the test facilities. The cost of testing will be considered to consist of the fol- lowing two components: Ci = The cost of repair of a unit which fails during test. When several independent failure modes are considered this definition can be expanded to include a different cost for each failure type. C2 — The cost per unit test time. 4. OPTIMIZATION BY DYNAMIC PROGRAMMING The profit from a unit test is the terminal value (given by the reliability value function) minus the test costs incurred during the testing process. At an intermediate time during the test the optimal test policy does not depend upon the costs already incurred since they are fixed or "sunk" costs. The optimal policy is that which maximizes the expected profit of the remaining portion of the unit test. The decision-making process is sequential; the decision to stop or continue is remade after each observation of the partial test results. Since partial test results are revealed continuously in time each time point is, in fact, a decision point. The basic repeated unit of the decision tree is shown in Figure 1. / / , DEC IS ION NODE STOP TEST FAILURE NODE -•-O CHANCE NODE NO FAILURE DECISION NODE FIGURE 1. The basic unit of the Decision Tree We define the optimal value function as V(6, t) =the expected value of the remaining portion of the unit test when the current state is (6, t) and when all future test decisions are optimal. Let us consider a small positive time interval, h, starting at a time when the state of the process is (6, t). The optimal value function, V(6, t), can be expressed in terms of the optimal decision at the start of the interval, the events which can occur during it, and the optimal value function at the possible resulting states. Let us consider the possibility of one event occurring in (t, t + A) at time r + e, where e < h. Since there are two possible decisions the optimal' value function is equal to the maximum of two terms. The equation is V(0, T) = Max{tf(0, t),-C 2 A + F (T + hlO)-Fo(Tldy 1-Fo(t/0) [-d + VlOniO, r+e), t h (6, r+e))] l-F (T+h/0) + l-Fo(rld) V(d,T+h) The term V(d, t) can be subtracted from both terms in the Max bracket. The resulting expression can be written as 0=Max (J(0, t)-V(0, t), h \-C 2 ,t),a{-< OPTIMAL BURN-IN TESTING 171 V(d, T+h)~V(d,T) + F (T+hld)~Fo(Tld) h(l-F o (Tl0)) [-C l + V(6 R {0,T + e),T R {d,T + e))-V(0,T+h)]\\. Let us now observe that we have the form O=Max{£/(0, t)-V(6, t), h£(h)}, where £(h) represents the interior of the second term. It should be noticed that the value of A has no effect through the term £(h) since when hi;{h) <0 the function V(6, t) is determined by the equation 0= U(0, t) — V(6, t). In fact, for any positive h, h£(h)<0 <^f(A)<0,Af(A)=0<^f(A)=0. Therefore, we can remove h and write the equation as O = Max{£/(0, t)-V{6, t), £(h)}. The interval, h, can be considered as an inspection interval with the understanding that if a failure occurs anytime during the interval, its occurrence will not be detected until time, t+ h, at which time the operation of the device will be resumed. If detection is not an important aspect of the problem, however, the limiting form of the equation is more useful. Provided only that the functions /o(t/0), 6h(6, t), and tr{6, t) are continuous functions of t everywhere in their range of definition, we can now let h approach 0. In the limit, we obtain the equation 0= Maxjf/(0, t) - V(6, r) , dV{ ^ T) - «o(t/0) [V(6, t) - V(d«, r«) + C,] -cX It is helpful to visualize the testing process as a point moving with time in the (0, t) space. As testing proceeds, r (effective age of the unit) increases and remains constant at non-failure points in time. When a failure occurs, both and r normally jump to a new value given by (0«, Tr). The solution of the dynamic program divides the (0, t) space into two types of regions: 1. Continuation region where (CI) U(6, t)-V(6, t) «0, and (C2) ^ ( ^ t) -/xo(t/0)[F(0, T)-r(0 ft ,T H ) + C,3-C 2 = O. 2. Termination region where (21) U(d, r)-V(d, t)=0, and (T2) d ^ T) - M o(T/0) [^(0, r) - F(0«, r«) + C,] - C 2 ^ 0. These regions are separated by a boundary, t*(0), which represents the decision to terminate testing. A simple illustration of a process trajectory in state space is given in Figure 2. 172 J. M. COZZOLINO STATE SPACE (ONE 9 PARAMETER) •(e) REPAIR JUMP TERMINATION BOUNDARY o Termination Point X Failure Occupied FIGURE 2. A process trajectory in state space The following theorem helps to characterize the boundary in terms of the properties of U(6, r). THEOREM: For each 6, the optimal policy, t*(0), is restricted to the following possibilities: a) t* = or oo. b) t* = t such that the derivative, 2 , is discontinuous at t. St , * . , dU(6, t) . c) t =t such that 2 is continuous at t, dr d U(0, t) 8V(6,t) , a 2 t/(6>, t) d 2 V{d, t) , and — s= — at r. 8t dr dr 2 3t 2 PROOF: Assume that, for given 6, V{6, t) and t*(^) are the optimal value function and decision boundary, and that a) 0<t*(6>) <oo b) U(6, t) has continuous first derivative, - 1 at t*. OT At a point {6, r* — e) (for e > 0) in the neighborhood to the left of t*, both V(d, t) and U{6, r) and their derivatives exist by assumption. Condition CI gives U(6,T*-e)-V(d,T*-e) ^0. Expansion about the point t = t* gives U(6, T*-e) - V(6, T*-e)=U(d, r*) - V(d, r*) -e dU(d,r*-ye) dV(6,T*-ye) dr St for some < y < 1. Conditions CI and T\ then imply the following condition in terms of the left-hand derivatives dV-{e,T*) ^ dU(e,T*) dr dr The left-hand derivative also satisfies Eq. C2 ay- (e, t*) dr fjLo(T*/e)[V(e, T *) - v(o R , t«) +c] -c 2 =o. In a neighborhood to the right of r* the conditions 71 and T2 give dU+(d, T*) dl ■ti (T*id) [U(e, T ) - v(o K , th) +C] -c 2 ^ 0. OPTIMAL BURN-IN TESTING Subtracting the latter two conditions and using 71, we have dV-(e,T) dU+(e, t) 173 dr dT Assumption C ensures that both one-sided derivatives of U are equal; hence we have the existence of Ps derivative at t* and the equality dV(e, r) dU(e, t) dT dr atr = r. Furthermore, if the functions U and V have second derivatives at t* then a Taylor series expansion using both derivatives and the above equation yields d 2 V(d, r) d 2 U(d, t) dT 2 dT 2 at r = r*. Hence part Q of the theorem is proved. The possibilities excluded in assumptions a and b are parts a and b of the theorem. The differential equation is a necessary condition upon the optimal value function. Its form and complexity will depend upon the nature of the 6 space and repair functions. In the simplest cases the complete family of solutions of the differential equation can be obtained analytically. The parameters V(8,r) for fixed 8, | 1 1 ___ L_- ---" V(8,t) = U(8,t-) V(8,r> _______ ; - ; ^^= , •' 1 _____ U(8,t) 1 | CONTINUATION REGION I TERMINATION REGION r*(8) SOLUTION AT TANGENCY V(8,r)FOR FIXED 8 CONTINUATION REGION | TERMINATION REGION 1 V = U ________ V(8,r> ^_____- :5 ==>-^' 1 ^^^Ui6,r) 1 1 1 SOLUTION AT DISCONTINUITY OF DERIVATIVE V(8,t)F0R FIXED 8 CONTINUATION REGION TERMINATION REGION V(8,r) ^ — - V=U U(8,r)-0 SOLUTION AT DISCONTINUITY OF U(8,r> Figure 3. Possible types of solution 174 J. M. COZZOLINO of its solution can then be evaluated using the above boundary conditions. This will be illustrated later for two special cases. Figure 3 illustrates the tangency condition and the two types of discontinuity conditions. The method of optimization discussed here is a continuous stage dynamic program. Other mathe- matically similar problems in the literature, such as that of age replacement policy for stochastically failing devices which wear out, have been formulated as a dynamic program using a "renewal equation" form. For example, [11], [15], [16], and [17] use that form. Either form could be used here although this differential form is thought to have some computational advantage for this problem. The Distribution of Remaining Test Time When analyzing the feasibility of burn-in testing it is useful to know certain statistics of the unit test. Perhaps the most important of these are the moments of the test duration which will be con- sidered here. Let t be the remaining test time when the test process is in state (0, t), D(t/6, t) its density, and D(S/6, t) be the Laplace transform of the density. An equation for D{Sld, t), similar to the value function equation, can be derived using similar argument. By considering a small time interval, h, and the events possible within it, the following equation is derived: D(Sld, t) = e-™ D(Sld H , t r ) ixti(rld)h + e- sh D(SI6, t) [1 -/lh>(t/0)A] for t<t*(9). By allowing h to approach zero, the limiting form obtained is dD(Sie,T) _^ {T/e) [D(Sldr) _ DiS i eRTH)]=SDiS/dT) for T<T * ( 0). The associated boundary condition is D(S/6,t*) = 1. As a consequence, the Arth moment of remaining test time, defined as i E k (d,T) = \ t k D(t/d, r)dt for it = 0, 1, 2, . . ., obeys the following equation: dE k (0, r) -/ao(t/0) [E k (d, t) -E k (0 H ,T R )] =-kE k .t(8, t) dT for A=l, 2, . . ., and t 1 t*(0). The boundary condition is E k {6, t*)=0 for A=l, 2, 3 The solution for these moments requires the knowledge of t*(0); then, since the (k+\) st equation involves the Ath moment, the equations for the moments can be solved successively starting with A=l. 5. EXTENSIONS OF THE FORMULATION The Optimal Policy Under Discounting In the analysis of processes which take place over an appreciably long time, thus tying up capital equipment and possibly delaying cash flows, it may be useful to account explicitly for "value of time" by optimizing the present discounted value of the process. Under this method, the values of all future OPTIMAL BURN-IN TESTING ' 175 value transactions are discounted from face value according to the length of time into the future at which they take place. Hence, if the amount, £/, is to be received at time, T, into the future, its present discounted value is Ue~®'\ where /3 is a positive discount factor. The optimal policy and value function under discounting will be different from that already dis- cussed; however, only a slight modification of the optimization equation is required to treat this case. The equation now appears as 0= Max {U(d, t) - V(6, T), dV( ° ,T) -pV(6, r) - u (tI6) [V(0, t) - V(0 H , r^-C,] + C 2 }- OT The one additional term is (3V(6, t). All that has been said about the non-discounted optimal solution applies to this case also. The effect of present value discounting will usually be to shorten the test. This can be expected since the terminal value is always received at the end of the test, whereas the test costs are spread out over the duration of the test. Therefore the expected present discounted value of a given policy is necessarily less than the expected undiscounted value of the same policy. Burn-In Testing Under Uncertainty of the Failure Process We have so far considered a known failure process. However, in many applications there exists a high degree of uncertainty about the underlying failure process. One method of relaxation of this assumption of complete knowledge of the density is to consider a known form of density function but with some unknown parameters. The state of knowledge about the values of such parameters is ex- pressed, at any given time, by a joint density function over all allowable sets of their values. This prior density is modified by the observation of results of the testing process according to Bayes' Rule. The concept of a conjugate prior density, as described by Raiffa and Schlaiffer [25], is used here to enable a mathematically tractable formulation of this problem. We suppose that the parameterization of the density includes parameter(s) cf> which are unknown, but do not change with failure. They are now included explicitly in the density, f(t/6, t, </)). The "prior" density, that which describes the state of knowledge at a given time, is denoted g((/)/i//) . It is normalized over the range of possible <f) values for any values of its parameters, t|/. The density of fife, unconditional upon the unknown 0, is Jail <* h(tie, r, <//) = f(t/e, r, 4>)g{<t>l<\>)d<i>. Jail <b The prior density is chosen to be conjugate to f(t/<f), t, </>) so that, after observation of a failure at time, t , following start-up of a unit of effective age, t, the posterior density, according to Bayes' Rule. /(tie, r, c/QgW) h(tld, t, i//) has the same functional form in as does g{<j), <//), but with parameters changed to «//f = t///. (6, t, (//, t). Hence the value of the parameters of the prior density is sufficient to describe the state of knowledge at any time. The failure rate function which appears in the equation is V(tld, T, (//)= /; h(x/6, t, ty)dx 176 J- M. COZZOLINO The state of the process under uncertainty is given by (6, t, i//, x) , where ^ andT are the values current at the time of last failure or start-up and x represents time since that event. Let this description be abbreviated by the definition d' = {6, t, ¥}. The state after failure in state (8' , x) is d' R {6\ x) = {6, t + x), Tr{6, t + x), V F {d, T, ¥, x)} x R (d', x) = 0. It is now apparent that the adaptive burn-in test problem, for one unit test, has the same form of opti- mization equation as previously described: 0=Maxj £/(0\ x)-V(0', x), dV{ ®' ,x) -v(xld')[V(0', x)-V(d' R , 0) -Ci]+&] . It will, however, be more difficult to solve because of the larger state space. 6. SOME SPECIAL CASES There are some important special cases of burn-in test processes for which there are no 6 type variables present; i.e., the failure density of a repaired unit is of the same form as that of some aged unit which never failed. In such cases, the effective age of the unit (t) completely describes the satte of the process and the differential equation holds on a line. We now consider three such cases which encompass a range of possible repair rule cases. Complete Replacement or Unrepairable Equipment (t« = 0) In the case of unrepairable equipment, a failure necessitates a complete replacement of tlv* .m J device by a new, untested one. Therefore the state of the process after failure is independent of the age attained by the old unit just prior to failure. This case is represented by the repair rule Tr(t) = 0. The resulting linear first order differential equation, ^^-- (jlo(tI6)[V(t) -V(0)]^ fil{Tl6)C 1 + C 2 is easily solved to give a family of all solutions parametric on one constant. For convenience, let w(t) =j- „ , -■ + C 2 f [l-Fo(xld)]dx Jo [l-Fo(Tld)] [l-F (Tld)] Since the constant of solution is determined by the boundary conditions we write it as t*. The solution is then ( V{t*) - W(t*) + W{r) for ^ r ^ t* V(t)=\ [U(t) fort*<r. The value of t* is then determined as indicated by the previous theorem. The independence of interfailure times allows easy solution for the Laplace transform of the density of remaining test time. It is D{S,t) = Q(s',t*) ''*'*'* f ° r T ^ T * OPTIMAL BURN-IN TESTING J 77 where 1-J e- Sx Mx/e)dx Q^,r)=— [ [_ Fo(Tle)] Although the solution for this case could easily be obtained by more direct arguments, it is pre- sented here both to illustrate the method and to enable a comparison between unrepairable and re- pairable devices. We can already see that burn-in testing is somewhat inefficient for unrepairable devices since a failure results in the loss of all test effort expended up to the time of failure. Superficial Failure (tr — t) In the previous case, the test time "invested" in the current unit was completely lost with a failure. It represents one extreme case the opposite of which is called here superficial failure and results in no loss of invested time. This case, represented by the repair rule t«(t) =t, is important because of the simplification of the decision structure. The differential equation for this case is The optimal value function in terms of t*, is \U{t*)-C, \ fXo(x/d)dx-C 2 (T-T) forr^r* V(t)=\ h [U(t) for t s*t*. The Laplace Transform for the density of remaining test time is easily found to be Z)(S/T)=e- s(T *-"forTS=T*. This is the transform of a unit impulse function located at t = t; consequently the length of remaining test time is not random. The sequential decision which remains unchanged in the face of partial test experience. At first glance it seems counterintuitive that a failure leaves the device no better or worse than if no failure occurred; however, there are reasonable assumptions which lead to this repair rule. Kooharian and Young [28] present a class of such models. Similarly, Cozzolino [7] discusses a special case of that class, based upon the concept of production defects initially present in the device. Bazovsky [4] gives a similar debugging model. Because of its special deterministic nature, the decision policy is easy to implement. Also, as we shall see, this model results in a relatively profitable burn-in test. Partial Renewal (tr = (xt). The repair rule tr(t) =oct for 0<a<l represents an intermediate degree of replacement. Consider, for example, a device composed of N indistinguishable and essential components. Failure results from failure of any one of the components. Repair is then accomplished by replacement of that failed component by a new, untested component. The test time previously accumulated on all other components is retained. This situation may be represented by a partial loss of effective age, where a in the repair rule is a fraction of the time retained; I J in this case. A more detailed discus- sion of this model is given in [7]. 178 J. M. COZZOLINO The differential equation for this class of repair rules is dV(i -£- L -VL (Tl8)[V(T)-V(aT)]=Vo(Tld)C 1 +C 2 . The combination of differential and difference operators present makes analytic solution difficult for most failure rate functions of interest. However, approximate solutions based upon numerical integra- tion are easily obtained. Since all solutions are within an additive constant of each other, the boundary conditions are easily satisfied. A Specific Example Some numerical results, based upon the above family of repair rules, were obtained in order to illustrate how the optimal test policy depends upon the repair rule parameter, a. The specific model used is based upon the density function, | b t a \ 6+1 Mtla, b)=- -3- *2=0, 6>0, a>0. The failure rate associated with this density is monotonically decreasing, b /xo(tla, b) = a + t The value of a tested unit is taken to be a function, R(t), of its first life after test, t, as shown in Figure 4. I 15 Ro' Ro' 85 Ro RO) ' VALUE RECEIVED FOR TESTED UNIT NORMAL INCENTIVE RANGE PENALTY RANGE RANGE Figure 4. Value vs performance The reliability value function is consequently the expected value of this function as of test termination and is a continuous function of age, U(t) -\: R(t)f(tl6,r)dt = R« 0.85 + 0.15 a + T a + T + 7 7 , + 0.15 a + r a + T + T 2 The parameter values chosen for the example are a = 72.0, 6=1.1, Ci = 1.0, C 2 = 0.01, R o = 90.0, ri = 500, and T 2 = 1000. The mean life is found to be 720 units and the value of an untested unit is U(0) = 78.57. The optimal test policy and value function were found for several values of repair rule parameter, a. The optimal policy, r*, is tabulated for these a values in Table 1, along with other quan- tities of interest. Figure 5 gives the optimal value function for one a value; a = 7/8. Figure 6 shows how the expected test profit, V(0) —U{0) increases with a. OPTIMAL BURN-IN TESTING 179 89 - 87 - ^ Sale Price 85 83 >Test Cost 81 79 / ■ >Tesl Profit 77 I 1 1 1 1 1 1 1 1 1 1 1 1 1 40 80 120 160 200 240 280 320 360 400 440 480 FIGURE 5. The optimal value function V{t) for a = 7/8 V(0)-U(0) 18 16 14 12 1.0 08 06 04 02 n - t — i — i— -~~\ i i i i i c ) 01 02 03 04 05 06 0.7 0.8 09 1 D 6 8 10 = l-a i iii mi i Figure 6. Test profit as a function of a; Case 3 The parameter, a, indexes a family of partial renewal models; the fraction of the device which is renewed with repair of a failure is (1 — a)= 1/iV. As a ranges from to 1 the model ranges from com- plete replacement (a=0) to the superficial failure case a— 1, and is therefore capable of representing a wide range of repair behavior. The computations show that as the degree of partial renewal decreases {a increases), the coefficient of variability of the total test duration decreases, as would be expected from previous results. As the test duration becomes more deterministic the optimal test policy, t*, lengthens. The expected profit due to testing increases and, in addition, the ratio of expected test profit to expected test cost increases, showing an inherently greater profitability as the renewed portion of the device decreases. This result conforms to our intuitive expectation based upon the previous results; as a increases a larger portion of the invested test time is retained at failure. Consequently, failure becomes less of a risk and the test process becomes closer to a deterministic one. 180 J. M. COZZOLINO Table 1. Optimal Test Policy for a Class of Repair Models a N=a-ar> Optimal policy T* Expected value Expected profit from test [V(0)-U(0)] Ratio (%) of test profit to test cost r V(0)-U(0) i IV(T*)-V{0)\ Expected duration of unit test Coefficient of variability of duration of unit test Before test V{0) After test V(T*) 1 14.5 78.595 78.979 0.0249 6.5 15.7 0.28 0.5 2 80.0 78.70 80.68 0.13 6.6 97.5 0.26 0.6 2.5 140.0 78.83 82.04 0.26 8.1 174.2 0.25 0.7 3.3 220.0 79.09 83.63 0.52 11.5 268.4 0.20 0.75 4 261.0 79.27 84.36 0.70 13.8 310.9 0.17 0.8 5 304.2 79.47 85.07 0.90 16.1 351.7 0.14 0.875 8 370.0 79.82 86.05 1.25 20.1 406.4 0.09 0.9 10 390.0 79.94 86.32 1.37 21.5 420.6 0.07 1.0 oc 470.3 80.43 87.36 1.86 26.9 470.3 0.0 7. CONCLUSIONS The analysis of the burn-in testing of repairable devices requires explicit information about the repair process in order to evaluate repaired devices. A class of repair models, having a form tractable with respect to the sequential decision problem for optimization of the burn-in test, is described by a repair rule of the form, 0«(0, t), tr(0, t), specifying the parameters of a density function within the family of the renormalized form of the new unit density. This family of densities has the same func- tional form of failure rate. The burn-in testing optimization, when treated as a continuous stage dynamic program, results in a differential equation similar to equations for the state probability of a stochastic process. In two important special cases, the differential equation is of the first order in total derivatives and conse- quently is easily solved and fit to the boundary conditions. In more complex cases it has a differential- difference nature. The wide differences in the expected test profit under different repair models of the partial renewal family leads to the conclusion that the explicit modeling of the failure-repair process is highly important for the burn-in testing of repairable equipment. While this result is based upon a specific class of models, it appears likely that the greater profitability of burn-in testing of repairable over unrepairable devices may prevail over a much wider class of models. REFERENCES [1] Barlow, R. E., A. W. Marshall and F. Proschan, "Properties of Probability Distributions with Monotone Hazard Rate," Annals of Math. Stat., 34, 375-389 (1963). [2] Barlow, R. E. and F. Proschan, Mathematical Theory of Reliability (Wiley, New York, 1965). OPTIMAL BURN-IN TESTING - 181 [3| Barlow, R. E. and E. M. Scheuer, "Reliability Growth During a Development Testing Program," Rand Memorandum RM-4317-NASA (Nov. 1964). [4] Bazovsky. L, Reliability Theory and Practice (Prentice Hall Inc., New Jersey, 1961). [5] Chemical Rubber Publishing Company, Standard Mathematical Tables, 12th edition. Commissioner's 1941 Standard Ordi- nary Mortality Table. [6] Cox, D. R.. Renewal Theory (Methuen and Company, Ltd.. London, 1962). [7] Cozzolino. J. M., "The Optimal Burn-In Testing of Repairable Equipment." Massachusetts Institute of Technology, Opera- tions Research Center Technical Report No. 23 (Oct. 1966). [8] Davis, D. J., "The Analysis of Some Failure Rate Data," J. Am. Stat. Assoc, 47, 113-150 (1952). [9] Drenick, R. E., "The Failure Law of Complex Equipment," J. Soc. Indust. Appl. Math., Vol. 8, No. 4 (Dec. 1960). [10] Dreyfus, S.. "Dynamic Programming and the Calculus of Variations," Rand Dept. R-441-PR. [11] Fox, B. L., "An Adaptive Age Replacement Policy." Operations Research Center, University of California, Berkeley 65-17 (PR) (1965). [12] General Electric Co., Transistor Manual, Chap. 18, "Semiconductor Reliability," 6th Ed. (1962). [13] Haddon, M. O, H. Cheesebrough, and R. E. Kirby, "Guaranteed Reliability at a Minimum Cost," IEEE Trans, on Reliability and Quality Control (Sept. 1958). [14] Hausman, W. H., and M. Kamins, "Early Failures in Automobile Parts: A Background Study in Reliability," Rand Mem. RM-4002-PR (1964). [15] Jorgenson, D. W. and J. J. McCall, "Optimum Scheduling of Replacement and Inspection," Operations Research II, 732-746 (Oct. 1963). [16] Jorgenson, D. W., and R. Radner, "Optimum Replacement and Inspection of Stochastically Failing Equipment," in Studies in Applied Probability and Management Science (Stanford University Press, Sanford, 1962), edited by K. J. Arrow, S. Karlin, and H. Scarf. [17] Kamins, M. and J. J. McCall, "Rules for Planned Replacement of Aircraft and Missile Parts," Rand Mem. RM-2810-PR (Oct. 1961). [18] Knight, W. R., "Exponential and Subsexponential Distributions in Statistical Life Testing," Ph. D. Thesis, University of Toronto (1958). [19] Lawrence, M. J., "An Investigation of the Burn-In Problem," Technometrics, 8, 61-71 (Feb. 1966). [20] Lomax, K. S., "Business Failures: Another Example of the Analysis of Failure Data," J. Am. Stat. Assoc, 49, 847-852 (1954). [21] Markowitz, D. H. and R. G. Pouliot, "Infant Mortality Screening; A Discussion of Burn-In and Other Screening Techniques," U.S. Naval Weapons Quality Assurance Office, Washington, D.C. (June 1967). [22] Morse, P. M., Queues, Inventory, and Maintenance (Wiley, New York, 1958). [23] Norris. R. H., "Run-In or 'Burn-In' of Electronic Parts," Proceedings of the Ninth National Symposium on Reliability and Quality Control (1963), pp. 335-357. [24] Proschan, F., "Theoretical Explanation of Observed Decreasing Failure Rate," Technometrics, Vol. 5, No. 3 (Aug. 1963). [25] Raiffa, H. and R. Schlaifer, Applied Statistical Decision Theory (Harvard University Press, Cambridge, 1961). [26] Riordan, J., Stochastic Service Systems (Wiley, New York, 1962). [27] Watson, G. S.. and W. T. Wells, "On the Possibility of Improving the Mean Useful Life of Items by Eliminating those with Short Lives," Technometrics, 3, 281-298 (May 1961). [28] Young. H. and A. Kooharian, "A Decision Model for the Operational Acceptance Test (OAT) of Minuteman Ground Elec- tronic Equipment," Sylvania Electronic Systems, Res. Rept. No. 477 (Sept. 1965). [29] Zelen. M. (Ed.). Statistical Theory of Reliability, University of Wisconsin Press, Madison, 1963. GOAL PROGRAMMING IN ECONOMETRICS W. Allen Spivey University of Michigan and Hirokuni Tamura University of Washington 1. INTRODUCTION Optimization models have been suggested as an alternative to simultaneous equation models for use in dealing with the problems of policy formulation and decision making at the level of the national economy (Theil [11]; Holt [5]; van Eijk and Sandee [14]; and Hickman [4]). In simultaneous equation models a policy maker prescribes a set of fixed values for goal or target variables and these are inter- related with a set of instrument variables by means of a system of simultaneous linear equations. The policy maker then solves this system (in which the number of equations or restraints is usually required to be equal to the number of unknowns) for the values of the instrument variables. It has been empha- sized (Tinbergen [12], [13]) that the policy maker's welfare function is not explicit in this model but is implicit in the selection of the goal and instrument variables. In the optimization approach one develops a welfare (preference) function which is optimized with respect to a set of econometric restraints. As is usual with such models, there can be many feasible solutions, and the optimization process then chooses a best or optimal solution from the set of feasible solutions. In the Theil optimization model (Theil [11]) a preference function is developed in which both goal and instrument variables appear. A desirable configuration of goal and instrument variables, not neces- sarily satisfying the constraints, is specified by the policy maker, and a quadratic function in deviations between the desired and feasible values of the variables is optimized (minimized). Although the quad- ratic function allows the policy maker to take account of a decreasing marginal rate of substitution between any two variables (at least over an appropriate subset of the domain of the function), it has the disadvantage of assigning the same penalty to overattainment of a goal as to underattainment (Hick- man [4], p. 6). A policy maker, for example, might not want to accept underattainment of a goal, but might be indifferent as to whether the goal was met or exceeded.! In this paper, we introduce an optimization model which differs from that of Theil in several ways. It permits overattainment and underattainment of any goal, yf, to be weighted either equally or dif- ferently in the decision maker's preference function. The model can be regarded as a special type of linear optimization model (as we make clear below) so it is easy to handle computationally and a wide range of parametric or sensitivity analyses can be implemented on a computer. These analyses can tStill another approach has been developed (Spivey and Tamura [9]), that of using a generalized simultaneous econometric model (GSEM), where the word generalized is used to emphasize that the number of equations or constraints is not required to be equal to the number of variables. The GSEM can be regarded as a link between conventional simultaneous equation models and the Theil model since, like the former it consists of simultaneous linear equations and like the latter a GSEM can be shown to have a solution which is identical with an optimal solution to a Theil type optimization model. 183 184 W. A. SPIVEY AND H. TAMURA display interactions between alternate optimal policies and goal attainment and can assist the policy maker in studying explicitly many of the trade offs between underattainment, attainment, and over- attainment of selected goals. Such sensitivity analyses are not available for either the Theil model or Tinbergen type simultaneous equation models. The model can also easily accommodate explicit restrictions on instrument variables in the form of nonnegativity and upper and lower bounds. The model is that of goal programming, originally developed by Charnes and Cooper in another context (see Charnes and Cooper [1]; also Charnes, Cooper, and Ijiri [2]; Ijiri [6]). We develop an appli- cation of this model using Klein's Model I of the United States economy as a basis (Klein [7]) and display in numerical form and interpret a sequence of decision making implications of several classes of the sensitivity analyses that can be performed with this model. 2. GOAL PROGRAMMING AND GOAL ATTAINMENT A reduced form of a linear econometric model can be written (1) y=Rx + s, where y is an n by 1 vector of goal or target variables, x is an m by 1 vector of instruments,/? is ann by m (real) matrix of impact multipliers, and 5 is an n by 1 vector of constant or additive terms. If m — n, this is the reduced form of a conventional system of simultaneous linear equations; if m ^ n, then (1) is an example of a generalized simultaneous econometric model, in which m = n is a special case (Spivey and Tamura [9]). Let y* be the goal vector prescribed by the policy maker and consider the model (2) y* = /?* + s, where R and 5 are assumed to be known. If R is singular, then (2) either has no solution or infinitely many solutions. If R is non-singular, then (2) always has a unique solution; however, it may not be acceptable to the policy maker because one or more components violate restrictions which are present in the mind of the policy maker, but are not explicitly a part of (2) — for example, nonnegativity.t We reformulate (2) as (3) y*-s = Rx-Iz + + h- x,z + ,z~mO, where / denotes an n by n identity matrix, R is the n by m matrix introduced in (2), and z + and z~ are n by 1 (unknown) deviation or discrepancy vectors. The ith components zf of z + and zf of z" measure the deviation upwards and downwards, respectively, of feasible y, values from the corresponding goal values yf; in other words, the ith constraint of (3) can be written (4) yt-s^RfX-zt + zi (»=1, ...,»), where /?, denotes the ith row vector of the matrix/? and s, the ith component of the vectors. Note that a nonnegativity restraint on x is explicit in this formulation. If the policy maker can find a feasible instrument vector x for which bothz+ andz,~ in (4) are zero, then the goal yf can be attained exactly. On the other hand, if there is no feasible vector x for which zt—zf = in (4) he cannot attain the prespecified goal yf. He can, however, find a feasible x that will tlf implicit restrictions on x are violated in the case R is nonsingular, Tinbergen [12], [13] would suggest model reformulation by selecting another goal vector y*' in an ad hoc manner, substituting y*' for y* in (2). and again solving for x. GOAL PROGRAMMING 185 allow him to come as close as possible to y* in a sense we will soon make clear. Now consider the following goal programming model (5) minimize e T z + + e T z~, subject to (6) y*-s = Rx-h + + Iz- x,z + ,z-^0, where e T is an n by 1 row vector t each of whose components is 1. The z+ and zf (i= 1, . . . , n) enter into the objective function (5) with a positive coefficient so the minimization process reduces the values of these variables in an optimal solution to the lowest feasible level and (5) attains a minimum value of zero if and only if both z + and z~ are the zero vectors (if and only if zt=Zi~ = Q for each i). The equations (6) clearly have feasible solutions because the rank of the matrix / assures that (6) can be solved for the z's alone. The advantage of reformulating (2) as a goal programming problem should now be clear: if (5) has a minimum value of zero, then (2) has a solution vector x which enables the policy maker to attain his prespecified goals yf (1 = 1, . . ., n) exactly. If it does not, then (2) is not solvable for a feasible instrument vector x. In the latter case, however, an optimal solution to (5) and (6) will disclose (a) how close the policy maker can get to the desired goal vector y*, and (b) whether this approximation is on the high side or low side for each individual goal, y* . The statement (b) can be made, and nonsense values of zf and zj avoided, because the vectors in the system of restraints (6) associated with zf and zj are negatives of each other. It is impossible for both zf and z,~ to assume positive values in an optimal solution because the basis property of linear independence is preserved in the simplex solution procedure. Denote an optimal solution to (5) and (6) as 'Xo (7) where xj= (*i,o v • ., x m ,o), Zq T = (zt, o> • • •'> z n, o) an( ^ z o T= ( z i~, o» • • •> z n , o)- There are exactly three possibilities in an optimal solution to (5) and (6): (i) zt,o = Zi,o — 0, (ii) <•> 0,^0=0, ( iH ) z<to=0, zf, >0. t The superscript T denotes a transpose. All vectors will be regarded as column vectors unless this superscript is used. The model (5) and (6) can also be written Minimize [0, e T , e r ] z* Subject to V y*-s= [/?,-/, /] lz< x z* \ 186 W. A. SPIVEY AND H. TAMURA Case (i) means that the ith goal y* is attained exactly whereas cases (ii) and (hi) indicate that a deviation, upwards and downwards respectively, must be accepted by the policy maker desiring a nonnegative solution to (2). The loose phrase, "how close the policy maker can get to the goal vector y*," is given a precise meaning in goal programming because in (5) and (6) one is minimizing W X \yf-(Rix+ Si )\, i where the vertical bars denote absolute value. An optimal solution x yields values of the discrepancy variables Z+ and zj which have the following effects: (i') Zi,o = z7,o =Q if y* = RiXo + sr, (ii') zt >0 =RiXo+Si — yf, z f " =0 if y?<R(Xo+Si; (iii') Zi,o=y*— (Rtxo + Si), zl o =0 i( y*> RiXo + si. This can be stated in still another way. Given an xo in an optimal solution (7), the corresponding target configuration y , which occurs when xo is used in (2), is (9) y = Rx + s. If case (i) occurs, yo = y* as desired, but in either case (ii) or (iii) y ¥= y*. In the latter cases yo is the closest vector to y* where closeness is measured by the absolute value metric or distance measure (this metric is sometimes called the L t metric). That is, yo is an optimal solution to the problem n (10) minimize^ \y*~ Yi\, subject to (11) y = Rx + s,x^0. The measure of closeness can be generalized by using an L p metric and minimizing (12) (2ly<*-*| p V subject to the restraints (11), where p is any integer greater than or equal to one. If p = 2, this is the same as minimizing the quadratic function X(y* — y;) 2 subject to the restraints (11), which is the form of the objective function in the Theil optimization model. The model in this case selects a feasible vector y which is closest to y* in the sense of Euclidean distance or ordinary least squares. The goal programming model (5) and (6) can be augmented so as to contain political and institu- tional as well as nonnegativity restrictions on the instrument variables x,(i = l, . . ., m). Such restrictions, called boundary conditions by Tinbergen ([12], [13]), must be dealt with implicitly in the use of his models, but can be incorporated as upper and lower bounds on xi. These offer no difficulty in a linear programming model. For example, we could (5) minimize e T z + 4- e T z~ subject to (6') y* - 5 = R X - Iz + + Iz- GOAL PROGRAMMING 187 Xi^ Li US, : k ^ U k keT, X, z + , z~ i=0, where S and T are any subsets of the set of integers {1, . . ., m} and Li and Uk are specified lower and upper bounds, respectively. The model (5) and (6') can be regarded as a decision model which explicitly takes account of all relevant restraints, both econometric and political or institutional. 3. WEIGHTING DISCREPANCY VARIABLES In addition to determining a feasible vector y which is closest to the desired vector y* in terms of the L\ metric, a policy maker might regard attainment of certain goals as more important than that of others and might wish to use different weighting factors on deviation variables to reflect this. One can accordingly replace e T in the objective function (3) by other nonnegative weighting vectors u T and v T . This gives the following objective function instead of (5), (13) u T z + + v T z~ which is minimized subject to (6). Also, u T = v T = e T is a special case in which each deviation variable is given a weight of one. Thus if a policy maker desires to attain a given goal yf more than others he can attach larger weights m and v ,• to the corre- sponding deviations zf and z,~, respectively. If he would accept an overattainment of yf, but cannot tolerate any underattainment, he could set Uj = 0, but would make Vj large enough to prevent zj from assuming a positive value in an optimal solution. t In this case, the model will allow the policy maker to impose the following condition (14) yf =i R jX + s h i.e., we can accommodate an inequality restraint if desired — a feature which is not incorporated in the Theil optimization model. Moreover, one can parameterize the objective function in a goal programming model to see how an initial optimal deviation package — say, one in which some goals are attained exactly and some are not attained initially — varies as different nonnegative weights are attached to the deviation variables. This can be contrasted with the effect of a weighting factor in a Theil model (L 2 metric model), which permits a delicate control of the size of the corresponding deviation symmetrically above or below the goal vector y* (for an empirical study of the difference in effects of the use of Li metric and L 2 metric, see Tamura [10]). 4. SOME EXAMPLES OF DECISIONS IN GOAL PROGRAMMING MODELS The classical Klein Model I (Klein [7]) will be used to present some numerical examples of goal programming decisions. Model I in reduced form is shown below; actual values have been given to the lagged variables and the trend variable so that model can be used in examining 1933 policy decisions. tt tFor any nonnegative vector of coefficients [u T v T ], (13) has a lower bound of zero in a minimization problem. Only models with nonnegative coefficients in the objective function will be considered in this paper. ttThis model was also used for illustrative purposes by Theil [11] and Spivey and Tamura [9]. 188 Consider (15) W. A. SPJVEY AND H. TAMURA 0.6659 -0.1883 0.6705' 0.2238 -1.2808 1.1188 0.1620 -0.2033 0.8102 0.0517 -0.2959 0.2584 0.6142 -1.4842 1.9289 0.0517 -0.2959 0.2584 + /37.58^ ' 9.0764 23.0109 -5.4938 32.0873 i 201.6060/ where C= consumption expenditures, P= corporate profits, W\ — private wages, /= private invest- ment, Y— national income, and K= end of year capital stock are endogeneous variables and W<l = gov- ernment wages, 7"= tax revenue, and G = government expenditures are exogeneous variables (all variables are measured in terms of billions of 1934 dollars). The constant terms are calculated under the assumption that the random disturbances are forecast accurately. EXAMPLE 1. Parameterization of Weights: Suppose that the goal vector selected by a policy maker is y T = (C, P, Wi, I) and the instrument vector is x T = {W2, T, G,). Then in the manner of Tin- bergen and Theil we can formulate a linear econometric model in these vectors y and x by extracting the first four equations from the reduced form (15). Let the goal values be set as a uniform 10 percent increase over the actual values of 1932 for each goal variable, i.e., y* T — (50.16, 7.7, 31.9, —5.58). Then we use the following goal programming model to examine 1933 policy decisions, (16) subject to (17) minimize zt + z c + zp~ + z p + z^ x + z Wl + zf + Zj + 0.6659X -0.2238 \ w -0.1620 I -0.0517/ -i\ , / 1 0/ 0.6705\ 1.1188 G + 0.8102 0.2584 U*i + W 2 , l , tr, Z c , Z c , Z p , Z p , Z Wl , Z Wi , Z/ , 2/ = U, where (16) and (17) correspond to (5) and (6), and where each goal is given equal weight. An optimal solution to this problem leads to the policy (18) 6.870 1 I 13.650|, ,15.770/ and the corresponding target configuration yo is (19) yo- GOAL PROGRAMMING 189 Comparison of yo with the prescribed vector y* indicates that the three goals C*, P*, and W* are attained exactly and the private investment goal /* is underattained by $.23 billion. Suppose that the policy maker, upon being informed of (19), wants to assign a larger weight to the attainment of the investment goal/* than is the case in the objective (16). Then the deviation variable zj can be given the weight A (z/ being given the same weight 1), and we (20) minimize z£ + z c + zp~ + z p + z+ 1 + z Wi + z/" + kz, subject to (17). We study this optimization for various values of k (continuous parameterization of k). Table 1 summarizes relevant information. For 1 Si k < 4.328 an optimal solution remains unchanged and for each such choice of weighting factor k, Xo in (18) is optimal and y in (19) is realized. This is indi- cated by the first data row of Table 1. For weight k ^ 4.328, we must shift to a new optimal solution x , getting a corresponding y . The second data row in Table 1 displays values of these vectors. The invest- ment goal is attained for weights k in the range indicated and we overattain the goal P* by one billion. Table 1. Optimal Policies and Consequent Goal Attainments Range of Values of A Components of y Components of xo Co Po r,, h w 2 , T Co lgK 4.328 50.16 LI 31.9 -5.71 6.870 13.650 15.770 4.328 S k < °° 50.16 8.7 31.9 -5.48 6.845 12.641 15.512 Note: Underlining indicates exact attainment of goal. Comparison of the two optimal policies indicates that attainment of /* is brought about by in- creased deficit spending; Go - To = 15.770— 13.650 = 2.120 billion for the original optimal x , whereas G — To = 15. 512— 12.641 = 2.871 for the alternate optimal xo- Note that the increased deficit is not attributable to increased values of Go, but rather to decreased values of T . That is, in this model optimal adjustments call for reducing taxes rather than increasing government expenditures. A policy maker might wish to inform himself further about the interactions between T, G, I, and optimal policies. T, G, and / are related by means of the fourth equation in (15), / = - 0.0517r 2 - 0.2959r+ 0.2584G - 5.4938. One could study how sensitive an optimal policy is to changes in either of the multipliers — 0.2959 and 0.2584. A computer program can be used to parameterize the matrix in (17), one row or one column at a time, so parameterizing any one element in the matrix is a special case. Specifically, one could determine a range of values containing — 0.2959 within which an xo remains optimal, assuming that other columns of the matrix are held constant. Similar remarks can be made concerning the para- meterization of the government expenditures multiplier 0.2584. Indeed, the following kinds of para- metric analysis can be performed on the model (16) and (17): (a) nonnegative parameterizing of the coefficients in (16); (b) parameterizing the column vector of constants in (17); (c) parameterizations of (a) and (b) simultaneously; (d) parameterizing the matrix in (17) one row or one column at a time. Thus the model provides the policy maker with considerable flexibility in both policy formulation and decision making. 190 W. A. SPIVEY AND H. TAMURA EXAMPLE 2: Prespecification of Both Goals and Policies: Theil [11] introduced an optimization model in which the policy maker specifies a goal vector, y*, and a desired instrument vector, x*. We reformulate this model by means of goal programming and introduce additional flexibilities by using alternate objective functions. The target or goal variables chosen by the policy maker are C, I, and D, where D is a distribution variable defined as W x — 2P (see Theil [11], p. 79). Prespecified goal values are, in billions of 1933 dollars, C* = 49.69, /* = — 3.10, and D* = 11.25. The desired instrument values for W 2 , T, and G are 5.04, 7.40, and 10.44, respectively (for specification of goal and instrument variables, see Theil [11]). Thus (21) y* T = (C*,I*,D*) = (49.69 -3.10 11.25) and (22) x* T = {Wt,T*,G*) = (5.04 7.40 10.44). The goal programming model is to (23) minimize e T z + + e T z~ subject to (24) /y*\ /s\ (R X, Z + , 2^0, , (A . where I in (24) comes from (21) and (22), and where / 37.5812 \ / 0.6659 -0.1833 0.6705' s = -5.4939 I, R= -0.0517 -0.2959 0.2584 \ 4.8580/ \ 0.2855 2.3583 - 1.4274, were calculated from data in Klein's model 1 in the manner suggested by Theil (see Theil [11] p. 87 and Tamura [10], p. 46), h and h denote identity matrices of size 3 and 6, respectively, and Z + = (Zf , Zj , z D , z W2 , Z T , Z (: ) , z~ = {zj}, ZJ, zjj, z W2 , Z T , z a ) , x T = (W 2 , T,G). As observed earlier, the objective function we are using assigns the same weights to each zj 1 " and zf. If the policy maker would accept an overattainment of goals C* and /* and an underattainment of goal D*, then minimization of the following objective function subject to (24) could be suggested, (25) f where/ is a column vector which contains 12 elements all but 3 of which are set equal to 1. The remain- ing three elements are the coefficients of zf , zf, and zj, and are set equal to zero. The zero coefficient of these variables then essentially excludes them from the minimization process. GOAL PROGRAMMING 191 Optimal policy solutions under the two different objectives are given in Table 2. For purposes of comparison we have also indicated a Theil type optimization decision, i.e., one from a L 2 metric model. TABLE 2. Alternate Optimal Policy Decisions (Billions of 1934 Dollars) Values of y Values of x C I D w 2 T G Desired Values 49.69 -3.10 11.25 5.04 7.40 10.44 Optimal Solution with (23) 46.35 -5.55 11.25 5.04 8.42 10.44 Optimal Solution with (25) 46.54 -5.25 8.84 5.04 7.40 10.44 Theil Model 47.77 -5.41 10.42 6.44 8.32 11.13 Note: Underlining indicates exact attainment of goal. 5. CONCLUDING REMARKS Earlier we indicated a generalization of optimization models in terms of the L p metric, a generali- zation which includes both goal programming and Theil type optimization models as special cases, i.e., subject to minimize I ^ \yf — J\ | p ) 1/p , y = Rx + s, x ^ 0. Although for p = 2 this is a nonlinear optimization problem and, as such, may appear tQ present insurmountable computational difficulties, recent research (Duffin, Peterson, and Zener [3] and Peterson and Ecker [8] ) has shown that this problem is one of a class of nonlinear optimization prob- lems known as geometric programming. Results for these models are available whereby numerical calculations and approximations can be developed. Moreover, certain types of duality theorems and "complementary slackness" properties have been established, the economic interpretations of which are yet to be explored. Research by the authors on these matters is now in progress. REFERENCES [1] Charnes, A. and W. W. Cooper, Management Models and Industrial Application of Linear Programming (John Wiley & Sons, Inc., New York, 1961), Vol. 1. [2] Charnes, A., W. W. Cooper, and Y. Ijiri, "Breakeven Budgeting and Programming to Goals," J. Accounting Research /, 10-43 (1963). [3] Duffin, R. J., E. L. Peterson, and C. Zener, Geometric Programming— Theory and Application (John Wiley & Sons, Inc., New York, 1967). [4] Hickman, B. G. (editor). Quantitative Planning of Economic Policy (The Brookings Institute, Washington, D.C., 1965), Chap. 1. [5] Holt, C. C, "Quantitative Decision Analysis and National Policy: How Can We Bridge the Gap?" in Quantitative Planning of Economic Policy, edited by B. G. Hickman (The Brookings Institute, Washington, D.C.. 1965). [6] Ijiri, Y., Management Goals and 'Accounting for Control (North Holland Publishing Company, Amsterdam; Rand McNally Company, Chicago, 1965). [7] Klein, L. R., Economic Fluctuations in the United States 1921-1941, Cowles Commission Monograph No. 11 (John Wiley & Sons, Inc., New York; Chapman and Hall, Ltd., London, 1950). 192 W. A. SPIVEY AND H. TAMURA [8] Peterson, E. L., and J. G. Ecker, Geometric Programming: Duality in Quadratic Programming and L v — Approximation, The University of Michigan, Technical Report No. 1968-2, 1968. [9] Spivey, W. A., and H. Tamura, "Generalized Simultaneous Equation Models." to be published in International Economic Review. [10] Tamura, H.. "Linear Models for Macroeconomic Policy Making" (unpublished Ph.D. thesis) The University of Michigan, 1967. [11] Theil, H., Optimal Decision Rules for Government and Industry (Rand McNally and Company, Chicago, 1964). [12] Tinbergen, J., "On the Theory of Economic Policy," of Contributions to Economic Analysis (North Holland Publishing Company, Amsterdam, 1952), 2nd ed., Vol. I. [13] Tinbergen, J., Economic Policy: Principles and Design of Contributions to Economic Analysis (North Holland Publishing Company, Amsterdam), Vol. XI. [14] van Eijk, C. J. and J. Sandee, "Quantitative Determination of an Optimum Economic Policy." Econometrica 27, 1-13 (1957). PRODUCTION AND EMPLOYMENT SCHEDULING IN MULTISTAGE PRODUCTION SYSTEMS Christoph Haehling von Lanzenauer* The University of Western Ontario London, Canada ABSTRACT A model is developed for planning optimal production and employment levels in multi- product, multistage production systems. The market requirements for each product over the planning period are given. Backorders and/or shortages are permitted. Backorders and shortages must be considered in order to determine the amount of each product's demand that should be filled, backlogged, or go unsatisfied if the production capacity is insufficient to fill all market requirements. Backorders and shortages, on the other hand, are desirable under certain dynamic market conditions. I. INTRODUCTION The production and employment scheduling problem — planning production and employment levels over time — has received a great deal of attention. Models and decision-rules generally based on single- stage production systems have been developed for many situations under various conditions [7]. When- ever many products are to be manufactured the concept of aggregate production scheduling has frequently been suggested [1, 5]. The aggregate production level is then represented by some appro- priate measure of production per unit of time [5] in view of the common denominator required. Once the aggregate production level is specified the problem of determining the quantities and the timing for the individual products still remains. This problem of disaggregation is of special importance when- ever the available capacity is insufficient to produce all market requirements and therefore some demand has to go unsatisfied. It is not at all obvious how the disaggregation decisions should be made in view of insufficient capacity. Following along the lines of more recent work in production planning [2, 3, 8], we consider a multi- stage production system consisting of various production facilities (stages) linked together by in-process inventories. The objective of this paper is to develop a model for multiproduct, multistage production systems by which optimal production decisions on a disaggregate basis for both sufficient and insufficient production capacity can be made. As a point of departure we take Zangwill's [8] model. The subsequent analysis, however, differs from Zangwill's work on a number of points: (a) The demand requirements need not be met. (b) Each production stage manufactures various different products. (c) Workforce decisions are integrated. (d) A different cost structure is assumed. II. THE MODEL The Production System The elements of a multistage production system are the individual production stages consisting of production facilities and an. inventory. These production facilities can be thought of as multipurpose * Assistant Professor of Business, School of Business Administration. 193 194 C. H. von LANZENAUER machinery. The output of a given production stage will be fed into the inventory bank from which the units can be used either as inputs to higher numbered production stages or to supply market require- ments for the stage's products. Figure 1 describes such a production system. "ff~ MARKET Production and Inventory Aspects Let Dj, be the market requirements of production stage fs (j= 1, 2, . . . , J) product i (i= 1, 2, ...,/) during the period t(t=l, 2, . . ., T). It is assumed that the requirements in each period over the planning period are known in advance. The planning period is medium range so no discounting is required. The market requirements are expressed in units. Note that the index i at each production stage represents physically different products. Let Hj, be the amount of product i in the inventory of stagey at the end of period t. The inventory equation can then be expressed by (1). (1) tih = H h + X ( +X j D h k=j+l n=l for all i, t andy= 1, ,y-i with XjtR {Xj l0 ) being the amount of product i produced in stage j during period t on regular time (overtime). aj£ (ajg 3= 0) can be considered as a production coefficient and is the number of units of production stager's product i required to produce one unit of stage &'s (k=j+l, j+2, . . ., J) product n (n= 1, 2, . . ., /). aj£ = O if product i of stagey is not being used for the production of stage Fs product n. The inventory equation for the last stage J of the production system is (2) Ui — JJi _l_ Vi _L Yi _ f)i for all i, t. The formulation of (1) and (2) assumes that the production capacity available is sufficient to satisfy all market requirements. If market requirements exceed the available production capacity, the excess demand must either be backlogged into the next period or has to go unsatisfied. The possibility for backlogs and shortages is only meaningful for market requirements since the requirements of sub- sequent production stages, i.e., the production levels, are variable themselves. The inventory equation, adjusted by backlogs and shortages, can then be expressed by (3) (3) n}, = Hj t _ l +X} at +X} lo -(Dl t + B} t _ l )+Bj t +Sl t - ^ z, -» *"*m > "tto k=j+] n=i for all i , t and y = (/' + 1 , J -I) with #j ( being the demand backlogged into period £+ 1 and Sj, the amount of the periods total demand (Dj t + Z?j,_,) which will never be filled. Note that the variables//],, Bj t , and SL are nonnegative vari- ables. Generally, one can assume that, whenever a positive inventory exists of production stage y's product i at the end of period t , Bj, and Sj, will be zero; on the other hand, whenever #j, and/or Sj, is greater than zero, there will normally be no inventory available. It should be realized, however, that the formulation of (3) permits the possibility of having positive inventory as well as backlogged and/or unfilled demand at the same time. This configuration can be desirable considering over time changing MULTISTAGE PRODUCTION SYSTEMS 195 economic conditions. For example, in view of a speculative increase in the selling price in periods ahead, it might be advantageous to let some of the current periods demands go unfilled, though inventory is available. The inventory equation for the last production stage is (4) //', = //),_, + X> m + X* M - (/>*,+ 0j,_, ) +B< ( + Sj, , for all i, t. For the special case of a production line (i.e., production stager's output will only be used in produc- tion stage;' + 1) we have instead of (3), (3a) ^ = ^- 1 +^-^o-(^ + ^- 1 )+^ + ^-2 «&, C*£„»+*W' n=l for aU i, t. In view of hmited in-process inventory capacities, we can restrict (5) EP/ffJ* 1 * for all;,*, i=\ where /3j is the amount of space required to store one unit of production stager's product i. Lj t repre- sents the available inventory capacity in space terms. Correspondingly, restrictions can be placed on the maximum amount to be backlogged per period. Such an upper limit is desirable for Bj T in view of the fact that the cost associated with backlogs are normally less than those associated with shortages. The production quantities are subject to the following capacity restrictions (6) X'W** 5 ^ for all;,, and (7) S 7 /*^ * for all;,,, where t) is the time required to produce one unit of production stage ;'s product i and Rjt (Ojt), the capacity in time units of regular time (overtime) available on stage; during period t. Workforce Considerations The production decisions can only be carried out if the necessary workforce is available. The workforce decisions are to be made in view of the production requirements and with respect to the dynamic costs involved. (For a discussion of dynamic costs see [5]). The workforce available during period t will be denoted by W t . According to [6] we can differentiate between "productive" (Wf) and "idle" workforce {W}). So This distinction is desirable in view of the possibility of keeping workers on a reduced pay rate rathei than laying them off. Let y) be the number of units of stage ;'s product i which can be produced bj one worker during period t. Assuming a constant worker productivity (i.e., a constant yj), the require< productive workforce for regular time or overtime production in t can be expressed by (9) ^-is^ii 196 C. H. von LANZENAUER and <io) ri-22$*ir The productive workforce in period t is therefore (id pf=£i^c^+*jw- The total workforce available in period £ can then be expressed by (12) £i^(^+^) + ^*i£^(^-«+^-io) + »7-i + ft-»-^-i for all*, with G?_i being the number of workers hired at the end of period t— 1 and Ff_i the number of workers laid off at the end of t— 1. Since the decision of which workers to lay off is usually made on grounds of seniority, we don't have to specifically account for this fact within the model. The same holds for the "idle" workforce which is generally determined by juniority. The number of workers to be assigned to the various production stages during period t can easily be obtained from (9) and (10). From an organizational point of view, it is necessary to maintain the proper relation between the number of processing workers and supervisory labor. Let Wf be the number of supervisory labor in period t, then we have an analogous situation to (12) (13) W t s =W t s . 1 + M l .i-N t - 1 for all f, with M t -\ being the number of supervisory labor added and N t -\ the number laid off at the end of period t— 1. The relation between the processing and supervisory labor can then be expressed as (i4) \t±j j V!m+xj»)'*rf. j=l i- with 8 being the span of control. The Objective Function The objective of the model is to determine the production and employment levels so that the con- tribution to profits over the planning period is maximized. Let r'j t be the selling price of one unit of production stager's product i during t. If we assume that r' is constant for all t , the total revenues of the production system are iiirjiDj-sj,), t = i j = i i=i since any amount backlogged into a future period will bring the same revenue. Assuming an over time changing selling price rj t , the total revenue depends on the actual price any backlogged demand can ultimately be sold for. If backlogged demand Bj t has to be filled at the existing price rk of period t the total revenues are t=i j=i i=i MULTISTAGE PRODUCTION SYSTEMS 197 If on the other hand backlogged demand will be filled at current prices, the total revenues become t=i }=i i=i Production costs are the sum of material and labor costs. Capital cost will be considered as fixed. The cost to produce stager's product i depends not only on the decisions made on stagey but also on the production and employment decisions made on the lower-numbered stages which provide input to stagey. However, if one handles the labor costs en bloc, the material costs of stagey 's product i are only the costs for additional raw material since the costs of the inputs from lower numbered stages are already accounted for. Let cj t be the cost of additional raw material to produce one unit of stage y's product i during period t. The total material costs are then t=\ j=i 1=1 The inventory costs are the product of the average amount in inventory and an inventory cost factor z). Therefore, we have j t=i j=i i=i ^ Backlog and shortage costs are t=\ j=i i=i with bf being the penalty for backlogging one demand unit of stage /s product i for one period, s! is is the expected cost in terms of loss of goodwill, over and above the lost revenue, of one unit of stage /s product i that cannot be filled. The costs associated with the labor-force are t =i j=i ,=i J with p being the regular-time pay rate of a processing worker per period, k the overtime premium, and A the pay reduction for the idle workforce. h(hs) is the hiring cost per worker (supervisor), f(fs) the cost of laying off one worker (supervisor), and ps the supervisor's pay rate per period. Combining the above expressions* leads to the objective function (15) / I t=\ j=\ i=\ (is) tp=z j 2 W-sp-^«+*y "Assuming an overtime constant price r' . 198 C. H. von LANZENAUER - hG, ~/F t - p s Wf - h s Gf -fsFf} -iii{«-(4+^-(4+*^. z=i j=\ i=i "■ N ry \ Yj ' -- z\{H\ +H i l )-b i .B i - (r?+sf)S* -p{\-\)Wl-h-Gt-f-F t -p s W s t -h s Mt-f s Nt}- Maximizing (15) with respect to (3) or (3a), (4), (5), (6), (7), (12), (13), (14), and the nonnegativity constraints for the variables is a problem that can be solved by linear programming.t III. CONCLUSION The purpose of model was to formulate the production and employment scheduling problem on a disaggregate basis in multiproduct, multistage production systems. The model determines the amount of the demand for each product that should be satisfied, be backlogged, or remain unfilled. These results are important whenever the production capacity is insufficient to produce all market requirements. On the other hand, the formulation can be very useful under certain dynamic market conditions. The model determines also how the productive workforce is to be assigned to the various production stages. The model seems to be applicable to integrated iron and steel manufacturing systems as well as to various chemical processes. The author wishes to thank Andrew Grindlay and Howard E. Thompson for valuable comments and suggestions. REFERENCES f 11 Buffa, E. S., Operations Management — Problems and Models (John Wiley & Sons, Inc., New York, 1968) [2] Elmaghraby S. E. and A. S. Ginsberg, "A Dynamic Model for the Optimal Loading of Linear Multi-Operations Shops," Management Technology, Vol. 4, No. 1 (1964). [3] Haehling von Lanzenauer, C, "Formulation of the Production Scheduling and Sequencing Problem by Bivalent Linear Programming," Management Science, forthcoming. [4] Hanssmann, F. and S. W. Hess, A Linear Programming Approach to Production and Employment Scheduling, Management Technology, Vol., No. 1 (1960). [5] Holt, C, F. Modigliani, J. Muth, and H. A. Simon, Planning Production, Inventories, and Work Force (Prentice-Hall, Englewood Cliffs, N.J., 1960). [6] Lippman, S. A., A. J. Rolfe, H. M. Wagner, and J. S. C. Yuan, "Optimal Production Scheduling and Employment Smoothing with Deterministic Demands," Management Science, Vol. 14, No. 3 (1967). [7] Silver, E. A., "A Tutorial on Production Smoothing and Workface Balancing," Operations Research, Vol. 15, No. 6 (1967). [8] Zangwill, W. I., "A Deterministic Multiproduct-Multifacility Production and Inventory Model," Operations Research, Vol. 14, No. 3 (1966). tNote that ^ X 2 ri P\t is a con stant. A MUTUAL PRIMAL-DUAL LINEAR PROGRAMMING ALGORITHM Milton Y. Harris Naval Research Laboratory ABSTRACT A new primal-dual linear programming algorithm is exhibited. A proof is given that optimal solutions to both primal and dual problems (when such solutions exist) are found in a finite number of steps by this algorithm. A numerical example is included to illustrate the method. 1.0 INTRODUCTION Since the publication of the primal-dual algorithm by Dantzig, Ford, and Fulkerson, in 1956 (see Refs. [3], [5], [7]), apparently little work has been done in this area (see Gomory and Balinski [1]). This paper presents a primal-dual algorithm which is similar to that of Dantzig, Ford, and Fulkerson in that it uses the complementary slackness concept and alternates between primal and dual problems; how- ever, it differs in that it is simpler, but requires a separate Phase I stage.* Essentially, this algorithm consists in generating and solving a sequence of subproblems which use the constraints of the original problem. The objective function at each stage, however, is determined by the solution of the problem of the previous stage. Whether this algorithm is any more efficient than, say, the revised simplex algorithm remains to be seen.t Before exhibiting the algorithm, we first establish some notational conventions and state some necessary results. 2.0 NOTATION AND BASIC THEOREMS For the sake of simplicity, and without any loss of essential generality, we will consider only problems of the form maximize z—c'x subject to Ax^b (2.1) x^O where A is an {nxm) matrix, b is an {nx\) vector with nonnegative entries, c is an (mxl) vector of criterion function coefficients, and x is an {mxl) vector of variables. The "primes" attached to symbols denote the transpose of the indicated matrix. All entries in /I, b, and c are taken to be given constants. We may state the dual to this problem as (2.2) minimize u—w' b *For an explanation of the two stages of the simplex method, see Dantzig [4], ch. 5. t A Fortran program has been written for the CDC 3800 computer to compare the efficiency of this algorithm with the sim- plex method, but as yet, no conclusive results have been obtained. 199 200 M. Y. HARRIS subject to w'A S= c' (2.3) w' ^0 where w is an (nxl) vector of variables. These two problems may be rewritten by introducing slack variables 5 and v for the primal and dual, respectively, and taking the transpose of 2.3 as follows: maximize z — c'x minimize u — w' b subject to subject to Ax + Is = b A'w — Iv = c (2.4) (2.5) {x,s)^0 {w,v)^0, where / is an appropriately dimensioned identity matrix in each case. We will digress here slightly to define formally the notion of complementary slackness, a concept which is central to the algorithm. DEFINITION: For given basic feasible solutions (x, s) to the primal and (w, v) to the dual, complementary slackness is satisfied if, and only if, i) w • s — and ii) x - v — 0. We conclude this section by stating four well-known results which are essential to our algorithm. Proofs are given by Dantzig in [4]. Suppose we have any two basic feasible solutions (x, s) for the primal (2.4) and (w, v) for the dual (2.5). Then LEMMA 1: w'b^c'x. LEMMA 2: (x, s) and (w, v) are optimal solutions for the primal and dual, respectively, if, and only if, w'b = c'x. THEOREM 1: (x, s) and (w, v) are optimal solutions for the primal and dual, respectively, if, and only if," complementary slackness is satisfied for these solutions. FUNDAMENTAL DUALITY THEOREM: If a feasible solution to 2.4 exists and z of 2.4 is bounded above, then optimal feasible solutions exist for both the primal (2.4) and the dual (2.5). 3.0 THE ALGORITHM The algorithm is initiated with basic feasible solutions to both the primal and dual, either or both of which may be contained in the intital formulation of the problem or may be generated by Phase I of the simplex algorithm. We also assume that z of 2.4 has a finite upper bound. The fundamental duality theorem assures us that optimal solutions for both the primal and dual exist. Recall the primal and dual constraints may be written as (3.1) Ax + Is = b (primal) and (3.2) w'A-v'I=c' (dual). MUTUAL PRIMAL-DUAL ALGORITHM . 201 Multiplying 3.1 by w' on the left and 3.2 by x on the right and applying Lemma 1 gives (3 3) ws + vx = wb — c-x^0. Let a = w • s + v • x. Now any decrease in the value of a brings us closer to the optimum since a decrease in a represents a movement toward the satisfaction of complementary slackness and a convergence of wb and c-x toward their common optimal value (see Lemma 2 and Theorem 1). We now proceed to define the algorithm. Suppose we have a basic feasible solution {x°, s°) to the primal and a basic feasible solution (w°, v°) to the dual with (x°, 5°) • {v°, w°)=x -v + s°-w > (i.e., the solutions are non-optimal). Let Y° be a current full-tableau matrix for the primal, F° be a cur- rent constant column for the primal, Z 1 be a current full-tableau matrix for the dual, and Z* be a current constant column for the dual. Consider the linear programming problem (3.4) minimize a — (x, s) • (v°, w°) subject to (3.5) Y°(x,s) = Y° (x, s) ^ 0. Suppose the optimal (if one exists) solution to 3.4 and 3.5 is (x 1 , s 1 ) and (zSs 1 )-^ ,^ )-^. Let Y 1 and YJ be the matrix and constant column, respectively, corresponding to this solution. If 6*^ = 0, then {x 1 , s 1 ) and (w°, v°) are optimal solutions to the primal and dual, respectively, by Theorem 1. If oto > (recall from 3.3 that a<f 3= 0), consider the problem (3.6) minimize a^ (x l , s 1 ) • (v, w) subject to (3.7) Z*{w i v)=Z\ (W, V) 52=0. Let the optimal tableau to this problem be represented by Z 2 and Z\ and the optimal solution be (w 1 , v 1 ); let (x^s 1 ) ■ (v 1 ,w 1 )=at As before, if af = 0, we are through. If not, we return to the primal and solve (3.8) minimize a 2 = (x, s) • (v 1 , w l ) subject to (3.9) ° (x,s) >0 for (x 2 , s 2 ), Y 2 , Y§, and a* = (x 2 , s 2 ) ■ (v 1 , w x ). The algorithm consists in following this procedure of using the optimal solution at each stage as the criterion function coefficients for the next stage until for some n, a* = 0. 202 M. Y. HARRIS In the next section we will show that the algorithm results in optimal solutions for both the primal and dual, whenever they exist, in a finite number of iterations. 4.0 PROOF OF CONVERGENCE We may state the subproblems described above in general as follows: Primal minimize ak = (x, s) ■ (v^ k \ w [k] ) subject to YM(x,s) = YM (4.1) (*, s)2=0, where [k] — k/2 for A; = 0, 2, 4, . . . and Dual minimize a m = (x lm \ s [m] ) ■ {v, w) subject to ZM(w, v)=ZW (4.2) (w, v) 5=0, where [m] = (m— l)/2+ 1 for m — 1, 3, 5, . . . We first show that the subproblems 4.1 and 4.2 have optimal solutions. THEOREM 2: Both problems (4.1 and 4.2 above) have optimal solutions for any ke{0, 2, 4, . . .} and any me{l, 3, 5, . . .}. PROOF: Both problems (4.1 and 4.2) have feasible solutions by assumption. Now for every basic feasible solution (x, s) to the primal and (w, v) to the dual, we have by Lemma 1 (4.3) x- v + w s^O. Since the problems in question are both minimizing problems whose objective functions are of the form of 4.3, we have just seen that both are bounded below by zero. By the fundamental duality theorem then, both have optimal solutions. Finally, we show that the algorithm does indeed converge. THEOREM 3: a„* = for some (finite) n. PROOF: We assume that the original problem has an optimal solution. Clearly a* < ao, since at worst (w°, v°) satisfies 3.7 and (x 1 , s 1 ) • (v°, w°) — c*o (this fact is just a general property of constrained optima, that is, that the optimal solution to a minimizing problem is no larger than any admissable solution). Suppose a* = a*. Then for every basic feasible solution (w, v) for the dual w • b — c • x x = w • s x + s l + x l • v (by 3.3) = «i? a* = ctQ (assumption) = w° • s 1 + x 1 ■ v° = w°-b-c-x 1 (by 3.3). MUTUAL PRIMAL-DUAL ALGORITHM 203 Thus for every basic feasible solution (w, v) for the dual w • b — c • x 1 2= w° • b — c • x 1 or w • b 2= w° • b, which implies that (w°, v°) is the optimal solution for the dual. Similarly, for every basic feasible solution (x, s) for the primal w° • b — c • x = w° • s + x • v° (by 3.3) = «i 2= a* — ct * (assumption) = w°-b-c-x 1 (by 3.3) Thus for every basic feasible solution (x, s) for the primal w° • b — c • x ^ w° • b — c ■ x l or c • x ^ c ■ x 1 which implies that (x l , s 1 ) is optimal for the primal. Now if a* > 0, then complementary slackness is not satisfied by the solutions (*', s 1 ) and (w°, v°) , and so at least one must be non-optimal by Theorem 1. This contradiction (i.e., we just showed that if a* — a*, both (m; , v°) and (x 1 , s 1 ) are optimal) im- plies that if a* > 0, then a* < a *. We may repeat this argument for ai and oc-z and so on. We conclude therefore that a * > implies that a* +1 < a* for each A; = 0, 1,2,. . .. Since there are only a finite number of basic feasible solutions to the primal and dual, and since each iteration of the algorithm results in a strict decrease in a, if optimal solutions exist, then by Theorem 1, there exists an n (finite) such that Furthermore the solution to the problem for which a n is the objective function is the optimal solution to either the primal or dual with the solution to the preceeding problem being the optimal solution to the other. 5.0 NUMERICAL EXAMPLE Consider the following primal problem (this problem was taken from Hadley [5], p. 269): maximize z= 2x x + 13* 2 + Tx 3 + 6x4 subject to 1 & 4*i + lx 2 + 3*3 + 5x 4 3 3= — xi + 3x2 + 3x 3 — 3x 4 (5.1) 7 3*2*1- 4x 2 +lx 3 + 8x4 5 5s 3xi + 5x 2 — 2x 3 — 4x 4 2 2= — 5xi — lx 2 + lx 3 + 6X4 8 5=6xi + 2x2-9x3+lx4 The dual is then minimize u = W\ + 3k;2 + 7w 3 + 5w 4 + 2w 5 + 8w 6 204 M. Y. HARRIS subject to (5.2) 2 =£ 4wi — 1^2 + 2tt>3 + 3u>4 — 5w$ 4- 6wa 13 =S 1 w\ + 3z#2 — 4m; 3 4- 5t^4 — \w$ + 2we 7 =S 3wi + 7w 2 + Iwa — 2w4 + lws — 9w 6 6 =£ 5^i — 3^3 + 8w 3 — 4w 4 + 6^5 + 1m;6 After appending slack variables S\ through se to 5.1 and v\ through y 4 to 5.2 and writing both systems in detached coefficient form, we have Primal (5.3) b Xi *2 x 3 *4 Sl «2 53 54 55 56 2 13 7 6 1 4 1 3 5 1 3 -1 3 7 -3 1 7 2 -4 1 8 1 5 3 5 -2 -4 1 2 -5 -1 1 6 1 8 6 2 -9 1 1 (5.4) Dual c W\ W-i w 3 Wi w 5 w$ fli V2 v 3 Va 1 3 7 5 2 8 2 4 -1 2 3 -5 6 -1 13 1 3 -4 5 -1 2 -1 7 3 7 1 -2 1 -9 -1 6 5 -3 8 -4 6 1 -1 Carrying out a phase / step of the simplex algorithm on 5.4 gives us the dual tableau (see Table 5.5) in basic feasible form (the objective row has been omitted since it is not needed in the rest of the procedure). Now we construct an objective row for 5.3 from the basic feasible solution to the dual (5.5). This is shown in Table 5.6. After adjoining 5.6 to 5.3, four pivots of the simplex algorithm are required to minimize. The result is shown in Table 5.7. Since the objective function has not been reduced to zero (entry (1, 1) in 5.7), at least one dual cycle is required. We construct an objective row for the dual (5.5) from the solution to the primal (5.7). This is shown in Table 5.8. After adjoining 5.8 to 5.5, one pivot of the simplex algorithm is required to minimize. The result is given in Table 5.9. Since the value of the objective function in 5.9 is now zero, 5.9 must be the optimal dual tableau (except for the objective row, of course). Thus the previous primal tableau (5.7) must have been optimal. Checking the values of the actual objective functions using the solutions given in 5.7 and 5.9 for the primal and dual, respectively, we find that both are 13. MUTUAL PRIMAL-DUAL /LGORITHM 205 ACKNOWLEDGEMENTS The author wishes to thank Dr. Richard Young, Associate Professor of Economics, Rice University, and Mr. Arthur Ziffer of the U.S. Naval Research Laboratory for their many helpful suggestions toward the development of this algorithm. Table 5.5 c w, W2 w 3 w 4 W 5 We V\ V2 Va v 4 1.2979 1.0000 0.0000 0.8511 0.0000 0.0000 0.3830 -0.0164 -0.0213 -0.0426 -0.0851 2.2098 0.0000 0.0000 -0.7328 1.0000 0.0000 1.0040 0.0247 -0.2044 0.0797 -0.0267 0.8179 0.0000 1.0000 -0.4206 0.0000 0.0000 -1.1544 0.0351 -0.0302 -0.1069 0.0421 1.8006 0.0000 0.0000 -0.0747 0.0000 1.0000 -0.0604 0.1227 -0.1336 0.0351 -0.0925 Table 5.6 b Xl x 2 X3 Xi 5l 52 5 3 Si 5 5 5 6 18.4018 0.0000 0.0000 0.0000 0.0000 1.2979 0.8179 0.0000 2.2098 1.8006 0.0000 Table 5.7 6 Xi X 2 x 3 X4 5l S2 S3 S4 $5 56 5.4018 0.0000 0.0000 -51.2537 0.0000 -3.4789 -6.2958 0.0000 1.8733 0.0000 0.0000 0.0000 1.0000 0.0000 -3.4930 0.0000 -0.0423 -0.4085 0.0000 0.2535 0.0000 0.0000 .0.0000 0.0000 0.0000 2.6338 1.0000 0.1972 0.2394 0.0000 -0.1831 0.0000 0.0000 11.0000 0.0000 0.0000 2.1268 0.0000 -0.7606 0.6479 1.0000 0.5634 0.0000 0.0000 1.0000 0.0000 1.0000 3.8028 0.0000 0.1831 0.4366 0.0000 -0.0986 0.0000 0.0000 3.0000 0.0000 0.0000 -28.4648 0.0000 -1.2113 -3.0423 0.0000 2.2676 1.0000 0.0000 6.0000 0.0000 0.0000 1.7183 0.0000 -0.3099 1.3380 0.0000 -1.1408 0.0000 1.0000 Table 5.8 c Wi U>2 W 3 Wa M>5 We V, V2 va V4 5.4018 0.0000 0.0000 11.0000 0.0000 3.0000 6.0000 0.0000 1.0000 0.0000 0.0000 Table 5.9 c W\ Ul2 w 3 W4 (#5 We Vl V2 Va V4 0.0000 0.0000 0.0000 -11.0000 0.0000 -3.0000 -6.0000 0.0000 -1.0000 0.0000 0.0000 2.8593 1.0000 0.0000 0.7863 0.0000 0.8672 0.3306 0.0000 -0.1372 -0.0122 -0.1653 1.8473 0.0000 0.0000 -0.7178 1.0000 -0.2013 1.0162 0.0000 -0.1775 0.0726 -0.0081 0.3028 0.0000 1.0000 -0.3992 0.0000 -0.2861 -1.1371 0.0000 0.0080 -0.1169 0.0686 14.6748 0.0000 0.0000 -0.6088 0.0000 8.1500 -0.4923 1.0000 -1.0888 0.2861 -0.7539 REFERENCES [1] Balinski. M. L. and R. E. Gomory, "A Mutual Primal-Dual Simplex Method." in Recent Advances in Mathematical Program- ming (R. L. Graves and P. Wolfe, editors) (McGraw-Hill Book Co., New York, 1963), pp. 17-26. [2] Dantzig, G. B., "Composite Simplex-Dual Simplex Algorithm-I," The RAND Corp., RM-1274 (Apr; 1954). [3] Dantzig, G. B., L. Ford, and D. Fulkerson, "A Primal Dual Algorithm for Linear Programs," Paper 7 in "Linear Inequalities and Related Systems," Annals of Mathematical Studies (Harold Kuhn and Albert Tucker, editors) (Princeton University Press, Princeton, 1956), No. 38, pp. 171-181. [4] Dantzig, G. B., Linear Programming and Extensions (Princeton University Press, Princeton, 1963), Ch. 5 and 6. 206 M - Y - HARRIS [5] Hadley, G., Linear Programming (Addison- Wesley, Reading, Mass., 1962), pp. 257-266. [6] Orchard-Hays, W., "A Composite-Simplex Algorithm-II," The RAND Corp., RM-1275 (May 1954). [7] Vajda, S., Mathematical Programming (Addison-Wesley, Reading, Mass., 1961), pp. 108-113. INTERCEPTION IN A NETWORK* R. D. Wollmer The RAND Corporation Santa Monica, California ABSTRACT This paper presents an algorithm for determining where to place intercepting units in order to maximize the probability of preventing an opposing force from proceeding from one particular node in an undirected network to another. The usual gaming assumptions are invoked; namely, the strategy for placing the units is known to the opponent and he will choose a path through the network which, based on this knowledge, maximizes his probability of successful traverse. As given quantities, the model requires a list of the arcs and nodes of the network, the number of intercepting units available to stop the opposing force, and the probabilities for stopping the opposition at the arcs and nodes as functions of the number of intercepting units placed there. From these quantities, the algorithm calculates the probabilities for placing the unit at the arcs and nodes when one intercepting unit is available, and the expected numbers of units to place at the arcs and nodes when multiple intercepting units are available. I. INTRODUCTIOIV The paper depicts a game between two opposing forces — an evader and an interceptor. Only the interceptor's problem will be dealt with. In this game, the evader attempts to proceed from one point to another in a network. The interceptor, who may possess one or more indivisible intercepting units, attempts to stop him by placing these units along arcs or nodes that he expects the evader to traverse. His problem is to place these so as to minimize the evader's probability of successful traverse. The evader's problem is, of course, to select a path that will maximize this probability. Both problems can be represented by a zero-sum two-person game. As the next section will show, however, the game matrix is very large and difficult to generate. By using a steepest ascent approach, a solution may be obtained much more easily for the interceptor. This paper presents an algorithm based on such an approach. As inputs, it requires a list of the arcs and nodes of the network, the number of intercepting units available to stop the opposing force, and the probabilities for stopping the opposition at the arcs and nodes as functions of the number (integer) of intercepting units placed there. From this information, the algorithm calculates the expected number of intercepting units to place at each arc and each node. For the case where only one inter- cepting unit is available, the expected values are probabilities of placement; for more than one intercept- ing unit, expectations are to be interpreted in the most obvious way. Specifically, if the expected number of units to place at a particular location is 4-1/4, one would always place four there and would place a fifth one there one-fourth of the time. The situation depicted by this algorithm is particularly relevant in applying interdiction to infiltra- tion and counterinsurgency, where the counterinsurgent or infiltrating force would play the role of the *This research is supported by the United States Air Force under Project RAND. Views or conclusions contained in this study should not be interpreted as representing the official opinion or policy of the United States Air Force. 207 208 R- D. WOLLMER evader, while the interdiction forces would play that of the interceptor. It also applies to police networks where a wanted fugitive (evader) is attempting to escape from a city while police (intercepting) units are attempting to close off his avenues of escape. The algorithm to be presented is a steepest ascent approach and yields solutions that are always optimal in the gaming sense (i.e., the evader's best chance of successful traverse is as small as possible) for the interceptor when he has but one intercepting unit at his disposal. For the case of multiple inter- cepting units, the solutions are not always optimal; thus for the general case the algorithm is an approximation scheme. II. INTERCEPTOR'S PROBLEM FOR ONE INTERCEPTING UNIT The network is characterized by sets of elements called arcs and nodes. Nodes are points or junctions, and arcs are undirected line segments joining nodes. An arc that joins node i to node j may be disignated by the symbol (i,j) or (J, i). The node from which the evader starts is called the source; the one he attempts to reach is called the sink. It is assumed that the evader and interceptor are each composed of one indivisible unit. Later on this restriction will be relaxed for the interceptor. The evader is attempting to choose a path from source to sink in such a way that his probability of successful traverse is maximized. The interceptor wishes to place his intercepting unit so that this probability is minimized. As is customary in game theory, it is assumed that each side's strategy will become known to the other and consequently each will need to adopt a mixed strategy. The inadequacy of pure strategies for the interceptor may be dem- onstrated by any network whose source and sink are relatively invulnerable to attack and any particular intermediate arc or node may be bypassed. Placing the intercepting unit at the source or sink would be relatively ineffective while placing it anywhere else would be completely ineffective as the evader would choose a path that bypassed the unit. Quantities defining the local effectiveness of placing the intercepting unit at arcs and nodes, and quantities expressing the interceptor's strategy are as follows: p(i) = probability that the evader will be stopped at node i, given that node i is on his chosen path and the interceptor places his intercepting unit there. p(i, j) = probability that the evader will be stopped at arc (i, j) given that arc (i, j) is on his chosen path and the interceptor places his intercepting unit there. 77 (i) — probability that the interceptor places his intercepting unit at node i. 7r(i,y) = probability that the interceptor places his intercepting unit at arc (i,j). The probability that the evader can successfully cross a particular node, t, is 1 — n(i)p(i) and that of crossing a particular arc (i,j) is 1 — Tr(i,j)p{iJ). Assuming these probabilities are independent* and the evader attempts to reach the sink by the path i,, . . . i n , his probability of successful traverse is n n-1 (1) ^ = E[ [ l-1T (ir)p(ir)] ][ [1— ir(ir, ir+l)p(ir, U+l)]- r=l r=\ This quantity, K, will be referred to as the value of the path ii, . . ., l». Given the p{i), p(i,j), tt(i), and v(i,j), the evader's problem is to find a source-sink path of maximum value. The interceptor's problem, given the p(i) and p(i, j), is to choose 7r(t) 2*0 and *This assumption is in fact not true. To illustrate this, suppose ir(i) = l/2 and n(j) = 1/2. Independence would require that the probability of the intercepting unit being at nodes i and j at the same time be 1/4 while in reality it is zero. However, the next section will show that this assumption leads to no inaccuracies for the one-intercepting-unit case. INTERCEPTION IN A NETWORK ._ 209 7r(i, j) 3 s such that the maximum value of all source-sink paths is minimized, subject to the con- straint %Tr(i) +2,ir(i,j) *£ 1. Note that if all p(i) and p(i,j) are strictly positive, any optimal strategy will be such that strict equality holds in this relationship. To see this, suppose the intercepting unit is not placed at any arc or node with probability p (i.e., 2/7r(t) + 2/7r(i, j) = 1 — p). Then increasing all ir(i) and Tr(i,j) by equal amounts totaling p yields an improved solution. III. A STEEPEST ASCENT APPROACH FOR ONE INTERCEPTING UNIT The problem of finding the n{i) and ir(i,j) can in theory be solved by linear programming. The linear program solution would include an optimal strategy for the evader as well as the interceptor. The value of the game would be equal to both the evader's maximum guaranteed probability of successful traverse, and one minus the interceptor's guaranteed probability of stopping the infiltrator. However, the game matrix for the program would require a column for every possible source-sink path and a row for each of the interceptor's pure strategies. The number of source-sink paths is exceed- ingly large even for relatively small networks and in later sections, when multiple intercepting units are allowed, the same may be said of the number of pure interceptor strategies. The algorithm of this section is a steepest ascent approach which avoids these enumeration problems; however, while it yields an optimal strategy for the interceptor and a value for the game, it does not yield an optimal strategy for the evader. Initially, all 7r(i) and ir{i,j) are assigned values of zero. Then they are increased by small amounts A7r(i) 2=0 and A-7r(i, j) 2=0. The proportions for the Air(i) and ^rr(i,j) are such that the additional intercepting unit placement probability, 2A7r(£) +XA7r(i, _/), divided by the decrease in maximum path value, is minimized as these two quantities tend toward zero. This is equivalent to maximizing the decrease in maximum path value per unit of additional intercepting unit probability allocation. The 7r(t) and ir{i, j) are then increased again in the same manner until they sum up to one. For a specific path, t'i, . . ., i n , increasing Tr(i r ) by A7r(i r ) replaces the factor [1 — ir(i r )p(i r )] in Eq. (1) by [1 — (rr(i r ) + A7r(£ r ))p(£ r )J and reduces the value of the path from K to K — &K, where 1 l-ir(»r)p(»r) J" Solving for the quotient of the additional intercepting unit probability allocation and the decrease in path value, the following expression is obtained: (6a) AK K ir(ir) P(lr) If intercepting unit placement probability were increased at an arc instead of a node, the same argument would yield an expression for the additional intercepting unit probability placement divided by the decrease in path value which is identical to the preceding one, with arc probabilities substituted for node probabilities. Specifically, for arc (i r , ir+i), this would be (3b) **&' ir+l) ~ 7T(i r , Ir+l) A£ Klp(i r ,i r+1 ) Note that in order to obtain any decrease in maximum path value, the A7r(£) and kir{i,j) must be strictly greater than zero along a subset, C, of arcs and nodes that interesects all paths of maximum 210 R- D. WOLLMER value. An efficient allocation would require that all members of C remain on maximum value paths, for otherwise, lessening the probability increase at an arc or node not on a maximum value path by a small amount would decrease SA7r(i) -+- 2A7r(i, j) , while not affecting the maximum path value. Then, of course, this small amount of probability could be redistributed among arcs and nodes on maximum value paths to get a strict improvement. Thus, the Air(j) and k-rr(i,j) assigned to the chosen C must be such that the A£'s are equal for each of its members. In other words, they must be proportional to 77(0 (4) p(i) b) , -7T(i,j). p(l,J) The constant — , where K is the maximum path value, was dropped in Expression (4), of course. The total increase in intercepting unit placement probability per unit decrease in the value of a maximum probability path would be: -n(i,j) p(ij) Thus, the problem is reduced to finding a subset of arcs and nodes intersecting all maximum prob- ability paths that minimizes Expression (5). Initially, all paths have value one since all ir(i) and 7r{i, j) are zero. Thus C is required to be a set blocking all source-sink paths. If the nodes and arcs are assigned capacities equal to the quantities of Expression (4a) and (4b), the problem of finding C is reduced to one of finding a minimum cut.* However, this is equal to the value of the maximum source-sink flow and can be found by the maximum flow algorithm of Ref. [3]. Note that once additional intercepting unit probability is allocated to the minimum cut, C, the Capacity of each arc and node in C decreases by its additional probability increase while the capacity of any other arc or node remains the same. Thus, the value of C will be reduced by SA7r(i) +SA7r(i, j) while all other cuts will decrease by amounts that do not exceed this. Hence, C remains minimum and the entire unit of probability may be allocated to C. The 7r(i) and ir{i, j) may therefore be solved for as follows: 1. Assign all nodes capacities of llp(i) and all arcs capacities of l/p(i, j) (since all 7r(i) and 7r(i, j) are zero). 2. Maximize flow from source to sink. 3. Let C be the minimum cut set and ^ its value. Set ir(i) =p(i)/V for ieC and Tr(i,j) =p(i,j)IV for (i,j)eC. Set all other 7r(i) and Tr(i,j) equal to zero. Note that since only one cut is obtained, each path of maximum value has positive intercepting unit probability allocated to only one of its arcs or nodes, and both its value and traverse probability are equal to the probability associated with that particular arc or node. Hence, the independence assumption does not lead to inaccuracies for the one-intercepting-unit case. Finally, the maximum path value is 1 — -p, where V is the value of C found in step 3. *Essentially, a cut set is a set of arcs and nodes blocking all source-sink paths. Its value is the sum of the capacities of its arcs and nodes, and the minimum value of all cuts is equal to the maximum value of the source-sink flow. See Ref. [3 J. INTERCEPTION IN A NETWORK 211 IV. MULTIPLE INTERCEPTING UNITS The steepest ascent algorithm of the last section may be modified to handle the case where there is more than one intercepting unit; however, the interceptor strategy obtained in some cases may not be optimal and the game value obtained may be approximate rather than exact. This will be discussed more fully later on. The quantities measuring the local effectiveness of placing intercepting units at arcs and nodes and those defining the interceptor's strategy must be redefined as follows: p(i)k — probability that the evader will be stopped at node i, given that node i is on his chosen path and the interceptor places k intercepting units there. pi.ii J )k — probability that the evader will be stopped at arc (i,j), given that arc (i,j) is on his chosen path and the interceptor places k intercepting units there. tt{i) = expected number of intercepting units the interceptor places at node i. ir(i,j) — expected number of intercepting units the interceptor places at arc (i,j), where p(i)o andp(i,j)o are both identically zero. Also let 7r(i')=7r(;) — [«■(*)] ir{i,j) =ir(i, ;') — [n(i, j) ] . Note that 7r(i) and 7r(i,j) are merely the fractional parts of 7r(i) and n(i,j). It will be assumed p(i)k and p(i,j)k are increasing in k and (6) a) p(i)k+i—p(i)k^p{i)k—p(i)k-i b) p(i, j)k+i—p(i, j)k < p(i, j)k~ p{i, j)k-u k > 1. The assumption of Expression (6) essentially states that if an additional intercepting unit is placed at any arc or node, it cannot be more effective in stopping the evader than any previous unit placed there. It is often satisfied in practice as the effects of overkill become more pronounced as the number of units at a location increase. The simplest example of this is the case where each unit placed at node i has independent probability p of stopping the evader. For this situation p(i)fc— 1 ~~ (1 — p) k and p(i)k+i~ p(i)k = p(l— p) k which is certainly decreasing in k, thus satisfy (6). The values of 7r(i) and n(i,j) of course are not sufficient in themselves to define the interceptor's strategy as any non-zero value of 7r(i) or ir(i,j) may often be realized in several ways. However, the assumption of Expression (6) assures that for a given value of -n ( i ) , the probability of stopping the evader at node i should he attempt to cross it is minimized by allocating [7r(t)] + 1 units at i with probability 7r(i), and [7r(i)] units with probability 1 — Tr(i). In other words, if 7r(i) = 3-1/2, tt{i) would best be realized by assigning three forces at i half the time and four half the time, as opposed to such a policy as allocating two forces half the time and five half the time. This is shown in the theorem to follow. First, define P{k, i} as the probability that k units are placed at node i. Note that 2,kP{k, i} = ir(i) and the probability of thwarting an evader's attempt to cross node i is 1 — 2,P{k, i}p(i)k. Minimizing this prob- ability is equivalent to maximizing 2,P{k, i}p(i)k. THEOREM: The quantity (t')2P{/:, i}p(i) k is maximized, subject to the constraint (ii)^,kP{k, i} = ir(i) when P{[ir(i)] + 1, i}=--n-(0 ,<P{[n(i)], i} = l-n(i), and all other P{k, i} = 0. PROOF: Let h = max {k/P{k, i}>0}, /= min {k/P{k, i} > 0}, and d=h-l. Let P{k, i} = P{k, i} be a set of values of smallest d which maximizes (i) subject to (ii). Suppose the theorem is false. Then h-l = d^2. Let m = min{P{h, i}, P{1, i}}. Set P{h, i} = P{h, i} - m and P{h-1, i} = P{h — l, i} + m. Then set P{1, i} = P{l, i} — m and increase P{/-f 1, i} by m. (If h — l = 2, then P{l+l} = P{l+l} + 2m.) Either P{h, i} or P{1, i} is now zero, decreasing d by at least one unit. 212 R- D. WOLLMER The increase in (/) is equal to m(— p(i)i + p(i)i+\ — p(i)ft + p(0/i-i)- Since (ii) is still satisfied, this is strictly negative from our choice of P{k, i}. However, since h > / + 1, P(i)h — p(i)h-i =S p(i)i+i—p(i)i - p(i)i + p(i)i+i— p(i)h + p(i)h-i 2= for a contradiction. Q.E.D. The same results, of course, hold for the Tr(i, j) of the arcs. The significance of this is that the search for an optimal interceptor strategy may be confined to strategies of the above type. When this is done, the 77(1) and Tr(i, j) are sufficient to specify the interceptor's strategy. As in Eq. (1) of Sec. II, if independence is assumed, the value of a path is still the product of the probabilities associated with its nodes and arcs; however, the probability associated with node i is now 7T(j) (1— p(i)[ir(l)]+l) + (1-17(0) (1— P(i)[n(l)]) or ^~p(i)[Mi)] — rr(i)(p(i)i7T(i)]+i—p(i)[Tr{i)]), and that with arc (i, j) is l^P(^y)[7r(i,j)] — 7r(l, j)(p(i, j)[7r(i,j)] + l— P(i, j)[n(i,j)]). If ti, . . ., in is a path of value K, and 7r(i r ) is increased by A7r(i r ), reducing the path value to K — &K, the expressions for K and K — &K, as before, differ only in the factor for node i r . Thus, for sufficiently small A7r(i r ), the expression for the quotient of the intercepting unit increase and path decrease is A7T(l r ) 1 f 1—P{ir)[n(i r )) 1 -TT{l r ) \ , AX K [p(ir)[n(i r )]+l— P(ir)[7r(i r )] and decreasing path value by increasing the number of intercepting units at an arc instead of a node yields the analogous expression, A7T(ir, ir+l) If 1 — P(ir, ir+l)[n(i r ,i r+l )) A£ If 1 — PKlr, lr + l)[ir(i r ,i r+l )] ... 1 = -p\ — — : : T. — : ^ ?rUr, ir+\) \ A I p{l r , lr+l)ln(i r ,i r + 1 )] + l—PKlr, lr+l) [n(i r , i r + 1 )] J As in Sec. Ill, intercepting units must be allocated along all arcs of some cut set, and the problem is reduced to finding a minimum cut set, which in turn reduces to one of finding a maximum source- sink flow. The node and arc capacities are now (7) a) C(l)=—7-^. 7TT 77 (I) P(l)[77(t)]+1-P(l)[7r(i)] b) C(l,J)= ,. ., /• -x TT\l,J)- Note that the c(i) and c{i, j) are no longer decreasing functions of 77(1) and 77(1,7) except over intervals whose endpoints possess the same integer part. Specifically, c(i) is a decreasing function of INTERCEPTION IN A NETWORK 213 7r(t) for < -n-(i) < 1, 1 s£ 7r(i) < 2, etc. However, c(i) is discontinuous and actually increases at all integer points for which (6a) holds with strict inequality. This follows since lim l—p(t) [„(!)] . ,., _ 1— ■p(Q[ir(i)]+i < 1— p(t)[7r(»)]+i ff(i) ^ 1 p(i)[n(i)]+l— P(i)[ir(i)] P(i)[n(i)]+1—P(i)[n(i)] P (0 [n(i)]+2 ~ P (0 [ir(i)]+l Similar results hold for c{i, j). Thus, when a minimum cut set is found and additional intercepting units are allocated along its arcs and nodes, that cut remains minimum provided the expected number of intercepting units on these arcs and nodes do not reach or exceed their next highest integer values. This requires that the steepest ascent approach be modified as follows. Starting with all 7r(t) and ir(i, j) equal to zero, capacities are assigned to the arcs and nodes of the network which are equal to the quantities defined by Expression (7). Flow is maximized from source to sink and the 7r(t) and Tr(i,j) corresponding to the arcs and nodes of the minimum cut set are increased by amounts A7r(i) and A7r(i, j) which are proportional to their capacities. The constant of proportionality, M, is the smallest possible such constant that will either increase some 7r(t) or Tr(i,j) to the next highest integer or will increase %Tr(i) + ?.TT{i,j) to n. If the latter does not happen, the new values of ir(i) and rr(i,j) are used to calculate new capacities and flow is maximized again to obtain a new cut. The process is repeated until %7r(i) +£77- (1,7) = n. Specifically, the algorithm for the general case of n intercepting units is as follows: 1. Set all 7r(i) = 0and all ir{i,j)=0. 2. Assign the nodes and arcs of the network capacities of c(i) and c(i,j) respectively, where c(i) and c(i, j) are as defined by Expression (7). 3. Maximize flow from source to sink and let C be the nodes and arcs of the minimum cut set found. A n Tit {1 _ tHO} 4. Compute M x = mm - 1 ... i€C c(i) 1VI2 — mm — — (i,j*C C(l,j) 5. Let M=min (M u M 2 , M 3 ). Set bn(i)=Mc(i) if ieC and &ir(i,j)=Mc(iJ) if (i,j)eC. Set all other A7r(t) and Air(i,j) equal to zero. 6. Set ir(i) = ir{i) + Att({) and ir(ij) = n(i,j) + bir(i,j) for all i and (i,j). If %ir(i) = £w<i,» < n, go back to step 2. Otherwise terminate. Step 1 is simply an initialization step. Steps 2 and 3 determine which nodes and arcs should have their expected number of intercepting units increased and the proportions of the increase. Specifically the expected number of intercepting units will be increased along the nodes and arcs of the minimum cut set found in proportion to their capacities. Steps 4 through 6 determine the proportionality constant M, makes the appropriate increases, and makes a termination test. Note that if M = M\, some -rr{i) will be increased by 1 — Tr(i), and since 77(1) is the fractional part of 7r(i), this will increase 7r(i) to an integer value. Similarly if M=?M 2 , some v(i, j) will be increased to an integer value. If M = M 3 , total expected intercepting units will be increased by n-XTr(i) —%TT(i,j) to n. In the latter case all n inter- cepting units have been allocated and the procedure terminates. 214 R. D. WOLLMER At termination the ir{i) and rr{i, j) represent the expected number of intercepting units to place at node i and arc (i, j), respectively. As mentioned earlier, the strategy will be to place [77(1)] + 1 forces at node i with probability ir(i) and [77(1)] forces with probability 1 — Tr(i). The maximum path value can be kept track of during the course of the algorithm. Specifically, if K is the maximum path value at the beginning of an iteration, then the maximum path value at the end of the iteration is K(l — M), where M is that found in Step 5. Of course, K= 1 at the beginning of the first iteration. This is easily proven by letting V be the value of C, the minimum cut set. Then substitute the generalized capacities (Expression 7a— b) for the specialized ones (Expression 4a— b) in Expressions (3) and (5). These replacements for Expressions (3) and (5) will have the same meaning in the multiple intercepting unit situation as Expressions (3) and (5) did for the single intercepting unit situation. The replacement for (5), which represents the total increase in intercepting unit place- ment probability per unit decrease in the value of a maximum probability path, is: ^Air(i) +SA-n-(t, /) = 1 y AK K However, 2A7r(i) +£A7r(i, j) =MV where M is as defined in Step 5 of the multiple intercepting unit algorithm. Thus MV = \y AK K or K-AK = K(l-M). Note that when there is more than one intercepting unit to be assigned, the strategy obtained is not necessarily optimal. To illustrate this, consider the network of Fig. 1. Figure 1. A simple network with node 1 the source and node 6 the sink. Suppose p(l, 2)i = 5/8, p(l, 2)2=1, p(2) fc = l — (1/4)*, and all other p(i)k and p(i, j) k equal 1 — (1/2)*. If there were two intercepting units available, placing both at arc (1, 2) with probability one will make the evaders successful probability of traverse equal to zero; however, the algorithm, as may be easily verified, will place both units at node 2 with probability one giving the evader a successful traverse probability of 1/16. For one unit, however, the algorithm will optimally place it at node 2 giving the evader a successful traverse probability of 1/4. N 4 8 FIGURE 2. Sample problem network with node 1 as source and node 9 as sink. INTERCEPTION IN A NETWORK 215 V. EXAMPLE For illustrative purposes, the algorithm is used to find the optimal placement of forces for the network of Fig. 2, with the results of each iteration. Ten intercepting units are available, and all arcs and nodes have identical probabilities associated with them. Thus, for each node, i, and each arc (i, /), we have: p( p( p( p( p( p( P( P( P( ) , = 0.6666 p(i,j)i = 0.6666 ) 2 = 0.8332 p(i,j) 2 = 0.8332 ) 3 = 0.8748 p(i,» 3 = 0.8748 ) 4 = 0.8852 p(i,y) 4 = 0.8852 ) 5 = 0.8878 p(i,y) 5 = 0.8878 )« = 0.8884 p(i,y) fi = 0.8884 ) 7 = 0.8886 p(i, j) 7 =0.8886 ) H = 0.8887 p(i, y)« = 0.8887 ) 8 = 0.8888 p(i,» » = 0.8888 p(0 ,0 = 0.8889 p(i,j) ,0 = 0.8889 The problem terminates in eight iterations. The results of these iterations follow: Iteration 1: All arcs and nodes have capacities of 1.5. A maximum flow is 1 unit along path 1,2,5,6,9. A minimum cut set is node 1 and its value is 1.5. Increasing intercepting units at node 1 to 1 gives a maximum path value of 0.33333. Iteration 2: The capacity of node 1 is increased to 2.0; all other capacities are unchanged from iteration 1. A maximal flow is 1.5 units along path 1,2,5,6,9. A minimum cut set is node 5 and its value is 1.5. Increasing intercepting units at node 5 to 1 gives a maximum path value of 0.11111. Iteration 3: The capacity of node 5 is increased to 2.0; all others are unchanged. A maximal flow is 1.5 units along path 1,2,5,6,9. A minimum cut set is node 9 and its value is 1.5. Increasing intercepting units at node 9 to 1 gives a maximum path value of 0.03704. Iteration 4: The capacity of node 9 is increased to 2.0; all others are unchanged. A maximal flow is 1.5 units along path 1,2,5,6,9 and 0.5 unit along path 1,3,5,8,7,9. A minimum cut set is node 1 and its value is 2.0. Increasing intercepting units at node 1 to 2 gives a maximum path value of 0.01853. Iteration 5: The capacity of node 1 is increased to 4.0; all others are unchanged. A maximum flow is 1.5 units along path 1,2,5,6,9 and 0.5 unit along path 1,3,5,8,7,9. A minimum cut set is node 5 and its value is 2.0. Increasing intercepting units at node 5 to 2 gives a maximum path value of 0.00927. Iteration 6: The capacity of node 5 is increased to 4.0; all others are unchanged. A maximal flow is 1.5 units along path 1,2,5,6,9 and 0.5 unit on path 1,3,5,8,7,9. A minimum cut set is node 9 and its value is 2.0. Increasing intercepting units at node 9 to 2 gives a maximum path value of 0.00464. Iteration 7: The capacity of node 9 increased to 4.0; all others are unchanged. A maximal flow is 1.5 units along path 1,2,5,6,9 and 1.5 units along path 1,3,5,8,7,9. A minimum cut set is arcs (5,6) and (5,8) and its value is 3.0. Increasing intercepting units to 1 at arcs (5,6) and (5,8) gives a maximum path value of 0.00155. Iteration 8: The capacities of arcs (5,6) and (5,8) are increased to 2.0; all others are unchanged. A maximal flow is 1.5 units along path 1,2,5,6,9 and 1.5 units along path 1,3,5,8,7,9. A minimum cut set is nodes 6 and 8 and its value is 3.0. Increasing intercepting units to 1 at nodes 6 and 8 gives a maximum path value of 0.00052. Algorithm terminates as all 10 units have been assigned. 216 R D - WOLLMER The final assignment is 2.0 intercepting units each at nodes 1,5, and 9; and 1.0 intercepting unit each at nodes 6 and 8 and arcs (5,6) and (5,8). REFERENCES [1] Durbin, E. P., "An Interdiction Model of Highway Transportation," The RAND Corporation, RM-4945-PR (May 1966). [2] Ford, L. R., Jr. and D. R. Fulkerson, "Flows In Networks," The RAND Corporation, R-375 (Dec. 1960). [3] Wollmer, R. D., "Maximizing Flow Through A Network With Node and Arc Capacities," Transportation Science 2, 213-232 (1968). [4] Wollmer, R. D., "Removing Arcs From A Network," Operations Research 12, 934-940 (1964). [5] Wollmer, R. D., "Stochastic Sensitivity Analysis of Maximum Flow and Shortest Route Networks," Management Science 9, 551-564 (1968). [6] Wollmer, R. D., "Algorithms For Targeting Strikes in a Lines of Communication (LOC) Network", Operations Research, 18, 497-515 (1970). THE FIXED CHARGE PROBLEM* David I. Steinberg Washington University St. Louis, Missouri ABSTRACT The fixed charge problem is a nonlinear programming problem of practical interest in business and industry. Yet, until now no computationally feasible exact method of solution for large problems had been developed. In this paper an exact algorithm is presented which is computationally feasible for large problems. The algorithm is based upon a branch and bound approach, with the additional feature that the amount of computer storage required remains constant throughout (for a problem of any given size). Also presented are three suboptimal heuristic algorithms which are of interest because, although they do not guarantee that the true optimal solution will be found, they usually yield very good solutions and are extremely rapid techniques. Computational results are described for several of the heuristic methods and for the branch and bound algorithm. 1. INTRODUCTION The fixed charge problem is perhaps one of the most intriguing problems of mathematical pro- gramming. Its structure is almost identical with that of a linear programming problem; yet, the one difference — the existence of the fixed charges in its objective function — has prevented the development of any extensive theory concerning its solution. There are two known exact methods of solution: total enumeration of the extreme points of the convex set of feasible solutions, and a mixed integer-continuous variable formulation. However, total enumeration is impractical for large problems, and, at this time, methods for solving mixed integer- continuous variable problems involving a large number of integer variables are not effective enough to be of much use, either. Therefore, heuristic methods for obtaining a "good" (although not necessarily optimal) solution have been developed. Some of these methods will be discussed in Section 2. Also presented is an exact algorithm, based on a branch and bound approach, which has been developed by the author, along with several new heuristic methods. Computational results are given in Section 4. 1.1 Formulation of the Fixed Charge Problem The fixed charge problem is a nonlinear programming problem of the form: (1) Minimize z = cx + dy subject to: Ax — b, (2) x ^ 0, This research was performed under AEC Research Contract No. AT(11-1)-1493. Submitted originally as Report No. COO-1493-15. 217 218 D. I. STEINBERG where A is an m X n matrix of real numbers, b is an m-vector in an m-dimensional space, E m , c and d are n-vectors in E", x is an rc-vector of unknowns, and y is an n-vector such that fO, if*, = ^ ^={i, if,;> - It is assumed that all dj ^ 0. The quantities dj are the fixed charges which are incurred if the corre- sponding Xj are (strictly) positive. If all dj = 0, then (1) and (2) reduce to a linear programming problem. It is known [7] that the optimal solution to the fixed charge problem occurs at an extreme point of the convex set of feasible solutions (2). Since the total number of extreme points is finite, the problem could be solved by evaluating the objective function (1) at every extreme point and taking the minimum of these quantities. However, as is the case with linear programming, the number of extreme points increases combinatorially with m and n. Hence, it is desirable to develop a method that requires that only a subset of the extreme points be examined in order to determine the optimal solution. Two such methods are the mixed integer-continuous variable integer programming formulation and the author's branch and bound algorithm. 1.2 A Mixed Integer-Continuous Variable Approach Hadley [7] presents the following formulation of the fixed charge problem and suggests that it be solved by use of Gomory's mixed integer-continuous variable algorithm. Minimize z = cx + dy, subject to: Ax = b Xj-Ujyj^O, 7=1, . . ., n 0<ft<l, 7-1 n yj integers x^O, where the Uj represent known upper bounds for the xj. If A ismXn, then the above formulation has m + 2n constraints. It should be noted, however, that the upper bound constraints on the yj need not be expressed explicitly, since no optimal solution would have yj > 1. Thus, the actual basis size required to use Gomory's algorithm is m + n. However, compu- tational experience with Gomory's algorithm has been very disappointing, even for some fairly small problems [7]; therefore, this approach is not widely used to solve fixed charge problems. 2. SOME HEURISTIC METHODS The fact that the optimal solution to the fixed charge problem lies at an extreme point of the convex set (2) has led to the development of some heuristic approaches which use adjacent extreme point techniques to obtain good solutions. All of the methods discussed below are basically adjacent extreme point methods. There are also several heuristic methods for solving the fixed charge transportation problem. These are presented in [2] and [3], but shall not be discussed here. FIXED CHARGE PROBLEM . 219 The following notations are offered prior to presenting the heuristic methods: Xb~ vector of basis variables, c B = vector of prices for variables in the basis, aj=yth column of A, B= basis matrix, yj — B~ l a,j, and From the theory of linear programming [6, 7] it is known that once a vector a,- is chosen to enter the basis, the vector a B r which leaves the basis is determined by XBr (4) — = Min { XBklykj , for y kj > 0} . Yrj k If all dk = (i.e., the problem is all linear) the greatest decrease in the objective function at each iteration is obtained if the vector to enter the basis, a,, is chosen so that X t (cj — Zj) =Min | — {c k — z k ), for y rk > [ k [yrk J When this expression is nonnegative the optimal solution has been reached. For the fixed charge problem (one or more dj ^ 0) criterion (5) must be modified to account for the fixed charges; the greatest decrease in z is obtained if a, is chosen so that (6) ej = Min {e fc , for y rk > 0}, k where (7) e k = — (ck~ z k ) + (d Br — d k ) Vrk Equation (7) must be modified slightly in the presence of degeneracy, in an obvious way. (It should be noted that the variable which leaves the basis, x Br , will, in general, be different for each,/). Because the objective function is concave, it is not true that when all e* 2s a global minimum has been reached. This stage may merely be a local minimum. That is, no adjacent extreme point yields an improved value of r, but it is still possible that some other extreme point of the convex set will yield a better value of 2. The methods discussed below differ in their approach to the problem of trying to obtain a better solution once a local minimum has been reached. The first three methods are those of the author. 2.1 Heuristic One 1. By using a Simplex Algorithm, vectors are removed from and entered into the basis according to (4) and (6), respectively, until all e,- 2= 0. Let the current solution be denoted by xo, with corresponding value of the objective function zo. Store xq. Define two parameters: ao = maximum number of consecutive iterations allowed to find a new extreme point with a value of 2 smaller than zo. /3o = maximum number of consecutive iterations allowed in which the objective function does not decrease. Also, define two "counters-," ai and /3i, which indicate the current number of consecutive itera- tions for which zo has not been improved, and in which the objective function does not decrease, respec- tively. Set ai = /3i = 0. 220 D. I. STEINBERG 2. Find the e, which yields the smallest increase in z and insert a, into the basis. Set ai = /3i = l. Proceed to Step 3. 3. Find the e, which yields the largest decrease in z and insert a, into the basis. Proceed to Step 4. 4. Compare the new z, z\, with zo: a) If Zi < Zo, return to Step 1. b) If z\ = Zo, but Xi 9^ Xo, set ai = ai + 1. If ai =£ ao return to Step 3. Otherwise terminate. c) If Z\ = Zo and *i = x , set «i = ai + 1 and /3i = /3i + 1. If ai =£ « and /3i =£ /3 go to Step 5. Otherwise terminate. d) If z\ > zo, set ai = ai + 1. If ai =£ ao return to Step 3. Otherwise terminate. 5. Perform /3i consecutive interations in each of which the vector which yields the largest increase in z is inserted into the basis. Return to Step 3. This algorithm attempts to find better solutions by moving away from the current local minimum in as few iterations as possible. The first such attempt (Step 2) actually increases the objective function as little as possible; if this attempt fails, the objective function is increased as much as possible (Step 5) with the hope that upon once again decreasing (Step 3) a better local minimum will be found. The efficiency of this method, as well as its accuracy, depends on the parameters ao and /3o. In general, ao should probably be somewhat larger than /3o. Since computational experience [4] has shown that the average number of iterations required to solve a linear programming problem with m constraints is about 3/2m, a reasonable range for a might be a e[m/2, m]. As the computational results of Section 4 indicate, much smaller values may also be used fairly successfully. Of course, economic considerations may also play a part in the determination of ao and j8o. 2.2 Heuristic Two 1. Same as in previous heuristic, except that the entire Simplex tableau corresponding to xo is saved. Define a as before, set /3 = n — m, ai = l, /3i=l, and proceed to Step 2 (assuming again that all ej 3= at xo). 2. Beginning with the tableau corresponding to xo, find the e, which yields the 03i)th smallest increase in z and insert a, into the basis. Set /3i = /3i + 1, and if /3i =S (J3o + 1) proceed to Step 3. Other- wise terminate. 3. Find the ej which yields the largest decrease in z and insert a, into the basis. Set ai = ai + l and if ai =£ ao, proceed to Step 4. If ai > a , proceed to Step 2. 4. Compare z\ with zo: a) If z\ > zo, return to Step 3. b) If z\ — zo, but jci t^ xo, return to Step 3. c) If Zi = zo and x\~ xo, set ai = 1 and return to Step 2. d) If z\ < zo, return to Step 1. 2.3 Heuristic Three Heuristic Three employs a Monte Carlo approach to find an extreme point with a better value of z than xo (where again at xo all ej 3= 0). In this method, extreme points within a certain "distance" from Xo are randomly examined. "Distance" here shall mean the number of basis changes required to move from one extreme point to another. Let 8 = maximum distance from xo allowed and let ao again be the maximum number of iterations allowed to find an improvement in z. The procedure is as follows: 1. Insert a.j into the basis, where ej is the maximum decrease in z. Continue in this manner until FIXED CHARGE PROBLEM 221 all e, S 5 0. Denote the current extreme point by x . Set a.\ =0 and proceed to Step 2. 2. Randomly select (according to the procedure described below) a vector within a distance 8 of Jto. Denote this vector by X\. Set a t = a t + 1 and proceed to Step 3. 3. Compare z\ with zo: a) If z\ > zo and «i 3= a , return to Step 2. If z x > zo and «i > a , terminate. b) If z\ ^ zo, return to Step 1, replacing xo and zo by *i and zi, respectively. In order to complete the description of the above algorithm, the procedure for the random selection of extreme points must be described. For simplicity, the vectors selected to enter the basis are chosen one at a time (and the corresponding vector to leave is determined by (4)). The basis transformation is then performed before the next vector to enter the basis is again chosen randomly. Since there are (n — m) nonbasis vectors, the random selection of a vector to enter the basis is performed by generating a random number from a uniform random distribution on the interval (0, 1), multiplying it by (n — m), and rounding the result to the next highest integer, say k. Then the A:th nonbasic vector is selected to enter the basis. The process is repeated 8 times, and the resulting extreme point is the first randomly generated extreme point within a distance 8 of Xo to be used in Step 2. An approach very similar to the above is proposed by D. Denzler [5]; the basic difference is that Denzler's method does not restrict its random search for an extreme point within a certain distance of the current best solution. Many iterations are thus wasted. Intuitively, it would seem that extreme points close to one which has a known good value of z should also yield fairly good values of z, and so it would appear to be a better procedure to restrict the search for better solutions to extreme points fairly close to Xo rather than attempt to search randomly the entire set of extreme points. 2.4 The Cooper-Drebes and Cooper-Olson Heuristics Cooper and Drebes [3] have developed two other adjacent extreme point heuristics, which differ from the above three methods mainly in that in the Cooper-Drebes methods the objective function is modified at certain stages in the algorithms; in particular, the linear costs cj,j= 1, . . ., n, are modi- fied according to (8) c j = c j + ^,x j >0 x j Cj=Cj, xj = 0, where "*/' implies the current value of variable Xj. At other stages in their heuristics, Cooper and Drebes remove from the basis the vector with the largest fixed charge (of the current basic variables) and insert into the basis the vector with the smallest fixed charge (of the current nonbasic variables). The two methods differ only in their decision rules for employing the various techniques discussed above. A modification of the Cooper-Drebes heuristics has been developed by Cooper and Olson [4]. The modification consists of randomly generating an extreme point a distance not more than 8 (where distance is defined as in Section 2.3) from the current local minimum, and then the Cooper-Drebes heuristic is again employed in an attempt to improve the solution. The random extreme point generation is called a "perturbation." The parameter 8 is allowed to vary from about m/3 to m/2. If the search following the last previous perturbation yields a local minimum worse than the best obtained so far, 8 is reduced; if the objective function value is the same as the best, 8 is increased; if a better value is obtained, 8 is not changed. 222 D. I. STEINBERG The computational experience reported on these methods [3, 4] indicates that they will yield the optimal solution a high percentage of the time and, moreover, can be expected to provide at least a near-optimal solution. 2.5 A Separable Programming Approximation Technique Another technique which yields a local minimum is to formulate the fixed charge problem as a separable programming problem [7]. Letting Uj represent known upper bound values of Xj, consider the following representation of the objective function (1): (9) Z = J£ [djXy +(dj + CjUj)k 2j ] , where d°) xj = \ 2J Uj and (11) \0j+klj+\2j=l and where no more than two \kj can be positive for a given j, and these can only be positive if they are adjacent. The complete formulation becomes: N Minimize z= ]£ [djkij+ (dj + CjUj)\ 2 j], subject to ^ aijUjkij — bu i = 1 , m AfcjSsO, k = 0, 1, 2 and j=l, . . ., n. When this problem is solved using a restricted basis entry Simplex method [7], the optimal solution will be a local minimum for the fixed charge problem. There is no evidence, however, that this procedure is any better than solving the fixed charge problem by using linear programming and ignoring the fixed charges. Moreover, no attempts are made to improve the local minimum obtained. 3. A BRANCH AND BOUND ALGORITHM In attempting to develop an exact solution method for the fixed charge problem, it would be natural to begin by investigating the possibility of developing an adjacent extreme point method. However, since the objective function is not always monotonic (that is, there may not always exist a monotonically decreasing adjacent extreme point path to the optimal solution) some other approach is called for. In this section a branch and bound approach is presented. (A general knowledge of the branch and bound philosophy will be assumed in the following discussion. The reader is referred to [1] for a general dis- cussion of branch and bound methods.) The algorithm finds the global minimum solution without having to examine all basic feasible solutions of (2). In the next two sections, the general approach will be presented and the computational details will be given in Section 3.3. FIXED CHARGE PROBLEM 223 3.1 The Tree of Solutions Suppose the following tree is generated: The nodes are numbered from left to right and from top to bottom, so that each horizontal "line" of nodes may be thought of as a level; that is, the first level consists of only node 1, the second level consists of nodes 2 and 3, the third level consists of nodes 4, 5, 6, and 7, etc. The constraints (2) will be required to be satisfied at every node. In addition, each node (except node 1) will be required to satisfy some additional constraints, of the form "xj=0" or "xj > 0." Since only basic feasible (extreme point) solutions are of interest, the latter constraints actually require variable Xj to appear in the basic solution, at a positive level. At each node (except node 1) one additional constraint is imposed, so that the solution obtained at any given node k must satisfy this new constraint plus the constraints imposed at every other node in the path from node 1 to node k. Thus, a node at level L must satisfy {L— 1) additional constraints. If a path is terminated when the additional constraints imposed by the node are such that: 1. There are m constraints of the form "xj > 0", or 2. There are (n — m) constraints of the form "xj = 0", or 3. No feasible solution exists when the current additional constraints are imposed, then it is clear that the collection of terminal nodes will contain all the basic feasible solutions to (2). That is, for terminal case 1, m ay's are required to be positive and hence a basis is uniquely determined. For case 2, (n — m) x/s must be zero, and hence the remaining xfs uniquely determine a basis; however, this basis may not yield a feasible solution, in which case terminal case 3 applies. Now, of course, the value of the objective function z can be computed for each extreme point and the minimal value obtained; however, it is desirable to reduce as much as possible the number of extreme points which must be examined. 3.2 The Algorithm The following notations will be introduced: S* — {j\ x j must be positive at node k}, Sj = {j\xj must be zero at node k}, ' z k — a lower bound of 2 at node k, and z* = the global optimal value of z. 224 D - ! STEINBERG The constraint set at node k thus becomes: / Ax = b (13) T k = x ^ Xj > 0, >eSf Xj = 0, jeS k x is an extreme point. Suppose that an upper bound estimate of z*, z, is obtained. Then a path can also be terminated at node k if (14) z k =s z, since the basic feasible solution corresponding to z has as good a solution as any of the solutions which could be obtained by continuing on the current path beyond node k. This procedure will cut down on the number of basic feasible solutions that need to be examined; its efficiency depends heavily on how good the bounds z k are and on how close z is to the true optimal value of z. Any heuristic method can be used to obtain an initial z\ thereafter, every time a path is terminated because a unique basic feasible solution has been obtained, the corresponding value of z, z, can be calculated and compared with the current z; if z is smaller than z, then z and its corresponding x replace z and its corresponding x as the current upper bound estimate of the optimal solution. A method for calculating z k is given in the next section. The branch and bound algorithm can be summarized as follows: 1. An initial upper bound, z, of the optimal value of the objective function is obtained by means of a heuristic. This bound and its corresponding solution vector x are stored. 2. The paths are then generated so that in branching from one node to another, the node at which the additional constraint is of the form "xj > 0" is selected as the next node in the path, until one of five path terminating rules occurs. The first three such rules are given on page 223. Rules 1 and 2 in- dicate that a unique solution has been reached (although it may not be feasible). Rule 3 pertains to infeasibility; that is, no basic feasible solution exists for the constraints (13). Rule 4 is the initial bound- ing procedure: A path is terminated if the lower bound z k of the objective function is as large as the current upper bound of z*; that is, if (14) holds, then the basic feasible solution corresponding to z has as good a solution as any of the solutions which could be obtained by continuing on the current path. Hence, the path is terminated. Rule 5 pertains to the fixed charges themselves. When the only fixed charges whose subscripts are not in either Si or S k are all zero, then as will be seen in the next section, the lower bound z k is the actual minimum value of the objective function (1) subject to the constraints (13). Hence, the path can be terminated (and if z k < z, then z k and its corresponding solution vector become the new upper bound). 3. When the last path has been terminated, the current lower bound z and its corresponding x will be the optimal solution. 3.3 Computation of the Lower Bounds In order to obtain the lower bound z k some initial computations are made. In particular, it will be useful to know the level of degeneracy that exists in (2); that is, how many basic variables can appear in a basic solution at a zero level. To obtain this information, the following n linear programming problems are solved: FIXED CHARGE PROBLEM . 225 Minimize Zj = Xj, subject to: x is an extreme point of the convex set (2) and Xj is in the basis. The method for solving these restricted linear programming problems will be presented later in this section. If the optimal solutions to (15) are denoted by xf, respectively, then if all xf > 0, no degeneracy exists; if any x*—0, then the number of such variables determines the maximum level of degeneracy which could occur. It shall be assumed that: (16) di^d-2^ . . . ^d„. In practice (16) is really no restriction, since the vectors may only have to be reordered so that (16) holds. Let S$ = {j\j4S* US*}, P—{j\xj can be in the basis at a zero level}, c?=pn (Sfusf), N k = the number of elements in S k : , and N k i = the number of elements in Q. Thus S k is the set of subscripts whose corresponding variables do not appear as an additional constraint at node k; the set Q consists of those subscripts whose corresponding variables could appear in the basis at a zero level at node k. Consider the following problem: (17) Minimize z\~cx subject to: xeT k . The desired optimal solution to problem (17), for any node k, lies at an extreme point of the constraints (2) and also satisfies the additional constraints imposed at node k. Denote this optimal solution by z k and let (18) 2§ = 2 dj. Further, let z| denote the value of the sum of the last (m — N k — N k i ) fixed charges. Define (19) z k = z k + z k + z k . In order for z k to be a legitimate lower bound for z*, it is sufficient to show that for any solution xeT k (i.e., any solution x which satisfies the constraints at node k), (20) cx + dy^z k . But, cx + dy=cx+ V g?j + Y djy-j, (21) ;eSf Ui since for./€Sf, yj= 1 and for^'eS^', yj— 0. By definition, ex 3 s z k . (22) 226 D. I. STEINBERG Since the maximum level of degeneracy possible at node k is N£, there must be at least (m — N£) fixed charges whose yj = 1 in any basic solution at node k; but,N k of these fixed charges are determined by S^ (and are included in z k ). Thus, there are at least (m — N k — Nq) remaining fixed charges whose yj must be unity. Since z k is defined to be the sum of the (tii — N* — Nq) smallest remaining fixed charges, it must be true that (23) 2 d M * 2 " its* Addition of the three relations (18), (22), and (23) yields the desired result (20). Hence z k as defined by Eq. (19) is a legitimate lower bound for the optimal solution which could be attained at node k. In order to complete the description of the algorithm it is only necessary to describe the procedure for solving the linear programming problems (15) and (17). Since (15) is merely a special case of (17) only the latter is considered below. Consider the following problem: (24) Minimize Z\ = ex (25) subject to: xeT k . The set of extreme points which satisfy the constraints (25) is a connected set; that is, an adjacent ex- treme point path exists between any two points in the set. Suppose an initial basic feasible solution to (25) is known, and that a Simplex algorithm is used to improve the objective function. Then, in a finite number of steps, a stage will be reached (call it Stage One) in which the only vectors which could be entered into the basis to yield a decrease in the objective function would require one of the vectors whose corresponding jeS k to leave the basis. (It is assumed throughout the discussion that none of the vectors whose corre- sponding jeS k are considered for entry into the basis.) Let the current solution at this stage be denoted by x, with corresponding i\ = ex. It is possible that some other extreme point £o exists such that zo < i\ , but that x and xo are not adjacent extreme points. In order to determine whether or not 5c is indeed an optimal solution to problem (17), the following constraint is appended to (25): (26) cx + x n + l = z l . The basis size is now increased to m + l, and the (m+ l)st basic variable for the current solution is x n + i = 0. Now, however, it is possible that one or more variables may be brought into the basis without removing any of the variables whose corresponding jeS k . If $ is not an optimal solution to (17), then a solution exists for which x n +i > 0. In order to determine whether such a solution does in fact exist, the following test for optimality is employed: Let R k denote the set of variables which are not in the basis for x (and whose j4S k ) . Then, for each xi in R k , first insert xi into the basis at a zero level, replacing x n + u and then enumerate the extreme points (of basis size m + 1 ) of the convex set determined by (25) and (26), by solving Maximize z = x,, subject to: Ax = b CX + Xn+l = Z\ (27) Xi — x„ +2 = xt x ^ FIXED CHARGE PROBLEM • 227 xj > 0, jeS$ Xj = 0, jeS k x is an extreme point, where x, is the current value of x,-, initially zero. At each step in solving (27) the minimum increase in Xi is taken. In this manner every extreme point of (25) and (27) which contains x,- will be examined. If at any stage x n +i > 0, the procedure can be terminated, and by beginning with this new extreme point, the objective function (24) can be reduced in the usual way until Stage One is again reached. If, on the other hand, the optimal solution to (27) for each x\ is zero, or no extreme points for which x n + \ > have been found after solving (27) for all XitR k , then x is the optimal solution to problem (17). This procedure will not, in general, be as lengthy as it might at first appear. If x is the unique optimal solution to problem (17), then the problems (27) will all have optimal values of zero, and one iteration per problem will indicate this. Moreover, if x is not in fact optimal, is should not require very many iterations to find a solution for which x n +i > 0. It should also be noted that the true minimal solution to problem (17) is required only when i\ > (z — z k — z k ), where z is the upper bound estimate, since if z.\ =£ z the path will not be terminated and decreasing z\ further is of no interest. Hence the above test for optimality for problem (17) need be applied at any node only when i\ > z. The question still remains: How is an initial basic feasible solution to (17) obtained? The procedure is as follows: Let W\ = {xj\jeS k , x } not in current basis} W<i = {xj\J€S%, Xj > in current basis} . Thus, those XjeWi satisfy the constraints imposed at node k, and those XjeWt do not; in order to remove those XjtWt from the basis, the following linear programming problems may be solved, for each xfiW 2 : Minimize z=Xj, (28) subject to: Ax= b x^O. When z = 0, either xj has been removed from the basis, or it is in the basis at a zero level and can be removed; if in an optimal solution z > 0, no feasible solution exists with Xj not in the basis. It should be noted that once a variable XjeW 2 has been removed from the basis, then XjeW\\ from the theory of linear programming it is known that if it is desired that some variables not appear in a basic solution, they may merely be excluded from consideration for entry into the basis; thus, pro- ceeding to remove the variables XjeW 2 from the basis while never considering the variables XjeW\ for entry into the basis will yield in a finite number of steps a basic solution in which W\ = S^ and W 2 is empty. Now, the variables XjtS'l must be inserted into the basis: Let Fi = {xj\jeS\ and xj > 0} Y 2 = { Xj \J€S k 1 and Xj = 0}. 228 D. I. STEINBERG For those variables XjeY 2 the following linear programming problems are solved: Maximize z — Xj subject to: Ax = b (29) x ^ XjeYi remain in the basis at a positive level x is an extreme point. Since the XjeY x are already in the basis at a positive level, the problems (29) can be solved by the method discussed earlier in this section. If the optimal solution to any of the above problems is zero, then no feasible solution to problem (29) exists with that variable in the basis at a positive level; if any solution exists for which z > 0, then variable xj has successfully entered the basis (it is not necessary actually to optimize the objective function). At this point, Xj will now be in Yi and the next variable XjeY 2 can be inserted in the basis via problem (29). In a finite number of steps, a stage will be reached at which Y\ = S\ and Y 2 is empty. If none of the variables XjeS% were considered for entry into the basis during the procedure of solving the problems (29), then an initial basic feasible solution to problem (17) has been obtained. It should be noted that at no stage of the branch and bound procedure will more than one variable not satisfy the additional constraints imposed by any node k. Thus, the amount of actual computation necessary to obtain the initial solution is much less than one might at first suspect. Moreover, the amount of computer storage required by the algorithm does not build up as the number of nodes examined increases, but remains constant throughout the procedure (for a problem of any given size). This advantage does not exist in most other branch and bound algorithms. 4. COMPUTATIONAL RESULTS The primary difficulty in finding fairly large problems to test heuristic methods is that the optimal solution must be known. Cooper and Drebes [3] have randomly generated 290 (5 X 10) problems with the following properties: (1) |ay|*£20, (2) l=£c jS £20, (3) 1 *£ dj *£ 999, (4) The density of A is 50 percent. The optimal solutions to these problems were obtained by total enumeration. In order to take advantage of these problems, larger problems were constructed by forming the A matrix as follows: ""/I, • • • A 2 ■ • • A = A k FIXED CHARGE PROBLEM 229 where the Ai are (5 X 10) matrices. The resulting problem is of size (5k X 10k) and its corresponding solution is known. Hence problems of varying size can be formulated to test the heuristic methods. Because of roundoff difficulties, some of the test problems are not included in the results. These are noted under the proper table. Since the heuristics tested are all adjacent extreme point methods, each iteration corresponds to a regular linear programming transformation. Thus, the primary means of comparison among different heuristic approaches is the number of iterations before the best obtained solution is found. Since these methods are not exact, the "best obtained solution" will not always be the optimal solution, and in some cases will be a function of the parameters related to the particular method. The approximate computer times (for an IBM 360/50) per iteration for the four groups of problems for both Heuristic One and Heuristic Two were (in seconds): 0.10, 0.18, 0.60, and 2.4. The Cooper and Drebes paper [3] reports computational experience with two heuristics: both heuristics were tested on a set of 253-(5 X 10) problems; 240 were solved exactly by the first method in an average of 75 total iterations per problem, and 245 problems were solved exactly by the second method in an average of 100 total iterations per problem. The three (5 X 10) problems which neither heuristic solved exactly were all less than 10 percent from optimality. No information on the number of iterations required to achieve the best solution is given. On a set of 70-(15 X 30) problems, the two methods were combined; the resulting heuristic solved 63 problems exactly with an average of 1200 iterations per problem. The nonoptimal problems were less than 3 percent from optimality. 4.1 Heuristic One Problems of size (5 X 10) ,( 15 X 30) , (30 X 60) , and (50 X 100) were used in the testing of Heu- ristic One. Of the 250-(5 X 10) problems, 235 were solved exactly with olq = 2 and /3o = 2. The average number of iterations required to achieve the final solution was 3 and the average total number of iterations per problem was 6.8. Table 1 contains a summary of the 15 nonoptimal problems. When ao was increased to 3 and /3o to 4, problems 2, 94, and 128 were solved exactly. No further increase in ao and /3o yielded any improvements in the 12 remaining suboptimal problems. Ninety (15 X 30) problems were attempted, and with a = 2, /3 — 2, the optimal solutions to 75 problems were obtained, with an average of about 5 iterations before the best solution was obtained and an average of 12 total iterations per problem. No further improvement resulted when ao was increased to 10; the average total number of iterations increased to 16. /3o was increased again to 15, with an improvement in only one of the suboptimal solutions-problem 43. The optimal solution is Z* = 8317.7 and the best previously obtained solution was Z = 9172.4; with /3 = 15, the solutions 8964.7, 8593.0, and 8525.8 were obtained. The average total number of iterations was 23. Table 2 contains a summary of the suboptimal (15 X 30) problems. 230 D. I. STEINBERG Table 1. Heuristic One-(5 X 70) Suboptimal Problems Problem Optimal Best obtained Percent number" solution solution from optimal 2 2830.4 "3317.6 14.7 57 4015.8 4017.9 0.05 70 3363.0 3468.3 3.1 86 2359.4 2380.7 0.9 93 2351.7 2436.6 3.6 94 2373.0 "2899.1 1.1 99 2418.0 2473.0 2.3 128 3459.0 » 3470.6 0.3 131 1543.0 1766.1 14.5 137 3593.1 3885.1 8.1 181 2826.2 3183.0 12.6 188 2417.1 2440.1 0.95 235 2080.9 2155.9 3.6 267 3072.2 3170.9 3.2 279 2323.4 2533.7 9.1 "Problems not attempted: 29, 36, 51, 55, 66, 67, 73, 75, 82, 88, 98, 117, 119, 120, 126, 130, 134, 141, 147, 158, 162, 164, 166, 176, 185, 201, 204, 206, 211, 215, 219, 226, 236, 237, 238, 248, 250, 260, 264, 266, 275. b Exact solutions obtained with a = 3 and /3o = 4. Table 2. Heuristic One-(15X30) Suboptimal Problems Problem Optimal Best obtained Percent number" solution solution from optimal 20 7676.4 7726.9 0.67 29 7522.6 7544.0 0.28 30 9452.9 9458.2 0.06 31 9201.7 9286.5 0.89 32 6666.4 7192.6 7.89 33 8738.3 8793.3 0.63 40 11,117.3 11,365.9 2.23 43 8317.7 8525.8 2.50 45 6785.6 6938.3 2.25 46 8517.7 8809.7 3.43 53 6497.6 6647.7 2.31 61 7651.0 8007.7 4.66 63 7836.0 7859.6 0.29 85 6307.2 6635.2 5.20 93 8592.7 9047.5 5.29 " Problems not attempted: 1, 44, 64, 68, 79, 89. FIXED CHARGE PROBLEM 231 Thirty-six (30x60) problems were attempted and 26 were solved optimally with a — 6 and /3 = 25. The average number of iterations before the best solution was obtained was about 11 iterations and the total number of iterations per problem averaged about 42. A summary of the suboptimal (30X60) prob- lems is given in Table 3. TABLE 3. Heuristic One-(30 X 60) Suboptimal Problems Problem Optimal Best obtained Percent number " solution solution from optimal 15 16,975 17,002 0.04 16 15,868 16,479 3.98 17 18,441 18,4% 0.30 23 15,303 15,748 2.91 25 16,886 16,987 0.50 27 13,960 14,110 1.07 31 14,944 15,301 2.39 32 17,350 17,440 0.52 37 14,716 15,181 3.16 43 13,946 14,274 2.35 ° Problems not attempted: 1, 5, 10, 20, 22, 24, 34, 40, 45. Nine-(50X100) problems were also tested, with a ~ 10 and /3o = 40; the results are given in Table 4. Table 4. Heuristic One-(50 X 100) Problems Problem number Optimal solution Best obtained solution Percent from optimal No. of iterations to best obtained solution Total No. of iterations 1 2 5 12 14 15 18 19 20 26,761 28,842 25,545 29,350 24,444 29,380 26,844 25,903 27,822 26,761 28,842 25,545 30,182 25,085 29,481 26,844 26,350 27,899 47 22 21 14 28 13 18 24 18 87 62 61 49 46 48 51 49 48 2.83 2.57 0.34 1.73 0.28 Average No. of Iterations 22.8 55.7 4.2 Heuristic Two The same sets of problems used to test Heuristic One were employed in the testing of Heuristic Two. In some cases, again, roundoff difficulties caused some problems to be discarded, and these problems are again noted in the tables. 232 D. I. STEINBERG Of the 268- (5X10) problems attempted by Heuristic Two, 255 were solved optimally. A summary of the suboptimal problems appears in Table 5. The average number of iterations required to achieve the final solution was 6, and the average total number of iterations was 15.2. TABLE 5. Heuristic Two-(5 X 10) Suboptimal Problems Problem Optimal Best obtained Percent number " solution solition from optimal 28 2993.4 3025.4 1.1 50 2383.9 2813.4 18.0 71 3058.3 3470.5 13.5 93 2351.7 2436.6 3.6 99 2418.0 2473.0 2.3 108 2414.8 2423.8 0.4 119 3467.6 3560.6 2.7 120 3579.1 3734.5 4.3 131 1543.0 1766.2 14.5 137 3593.1 3885.1 8.1 154 2589.5 2622.7 1.3 245 3574.7 3635.8 1.7 279 2323.4 2533.7 9.1 "Problems not attempted: 23, 27, 30, 52, 61, 67, 76, 81, 94, 98, 110, 128, 149, 158, 165, 202, 210, 224, 225, 239, 249, 254. Of the 287- (5X10) problems attempted by at least one of Heuristic One and Heuristic Two, 280 were solved optimally; the following 7 problems were not solved optimally by either method: 93, 99, 119, 120, 131, 137, and 279. Eighty-four (15X30) problems were attempted by Heuristic Two, and 74 were solved optimally. The average number of iterations required to obtain the final solution was 12, and the average total number of iterations per problem was 43. Table 6. Heuristic Two-(15 X 30) Suboptimal Problems Problem Optimal Best obtained Percent number ° solution solution from optimal 1 8096.2 8142.6 0.57 31 9201.7 9286.5 0.89 36 6162.5 6171.6 0.15 40 11,117.5 11,365.9 2.23 43 8525.8 8537.4 2.50 44 6807.3 6905.5 1.44 46 8517.7 8809.7 3.43 52 7880.0 8214.1 2.97 79 7765.0 7840.0 2.45 93 8592.7 8803.0 5.29 "Problems not attempted: 16, 18, 21, 22, 26, 27, 33, 37, 61, 72, 80, 83. Of the 96 (15 X 30) problems attempted by at least one of Heuristic One and Heuristic Two, only the following 5 problems were not solved optimally: 31, 40, 43, 46, and 93. ,u i^riAKGE PROBLEM 233 TABLE 7. Heuristic Two-(30 X 60) Suboptimal Problems Problem Optimal Best obtained Percent number solution solution from optimal 1 15,490 15,536 0.30 20 18,407 18,655 1.34 22 15,333 15,535 1.32 23 15,303 15,595 1.91 26 16,848 17,346 2.91 29 15,773 15,889 0.73 34 16,912 17,076 0.97 41 15,601 15,662 0.31 45 14,243 14,267 0.17 47 15,852 16,062 1.32 i "Problems not attempted: 5, 9, 11, 13, 14, 16, 17, 19, 31, 36, 40, 42, 43. Thirty-six (30 X 60) problems were attempted by Heuristic Two, and 26 were solved optimally. The average number of iterations required to obtain the final solution was 25, and the average total number of iterations per problem was 81. Of the 46 (30 X 60) problems attempted by at least one of Heuristic One and Heuristic Two, 42 were solved optimally; the following problems were not solved optimally by either method: 1, 23, 34, and 45. Heuristic Two was not tested on any larger problems. 4.3 The Branch and Bound Method Fifteen (5 X 10) problems and 10 (15 X 30) problems were solved by the branch and bound pro- cedure. The results are given in Tables 8 and 9. The total possible number of nodes that might have to be Table 8. The Branch and Bound Algorithm-(5X 10) Problems Problem Number Time number of nodes (minutes) 1 66 0.23 2 24 0.16 3 20 0.16 4 44 0.26 5 20 0.13 6 14 0.11 7 20 0.13 8 74 0.31 9 10 0.09 10 38 0.23 11 34 0.17 12 12 0.10 13 40 0.16 • 14 24 0.13 15 42 0.30 Average 32 0.16 234 D. I. STEINBERG TABLE 9. The Branch and Bound Algorithm-(l 5X30) Problems Problem Number Time number of nodes (minutes) 6 1242 17.9 8 672 9.1 9 3144 47.3 14 1134 18.3 18 308 8.5 29 2294 43.0 34 502 26.4 35 334 4.8 37 904 11.6 38 1546 24.4 Average 1208 21.1 examined is 2" +1 — 1. The average number of nodes that were examined was 32 and 1208, for the two groups of problems, respectively. Also included in the tables is the computer time required to move through the nodes; these times are for an IBM 7072. It is felt that these times could be greatly improved if the computer code were made more efficient. There are many improvements that could be made in the code, such as rewriting it in machine language and utilizing a more efficient Simplex method (such as revised Simplex). The average number of nodes examined for the (5 X 10) problems was 2 5 out of a possible total of approximately 2 10 nodes; for the (15 X 30) problems the average number of nodes actually examined was approximately 2 10 out of a possible total of approximately 2 30 . Hence, the ratio of the exponents of the number of nodes examined to the number of total possible nodes decreased when the number of constraints was tripled (from 5 to 15). If it is assumed that this ratio continues to decrease linearly, then the number of nodes that would have to be examined for a problem of size (45 X 90) is about 2 20 nodes. Hence, it appears that the number of nodes examined grows rapidly with the number of variables, and so the efficiency of the computer code is extremely important. 5. SUMMARY AND CONCLUSIONS The computational results of Section 4 clearly indicate that the heuristic algorithms are orders of magnitude faster than the branch and bound algorithm. Moreover, the heuristics almost always provide very good, if not optimal, solutions. When two or more heuristics are combined, the likelihood of finding the optimal solution increases considerably. In fact, when the two heuristics of Cooper and Drebes [3] are combined with Heuristic One and Heuristic Two, all 290 of the (5 X 10) problems were solved optimally, only one problem of the 96 (15 X 30) problems (problem 53) was not optimally solved. The branch and bound algorithm presented herein is a new, exact algorithm for the solution of the fixed charge problem and appears to be much better than the mixed integer-continuous variable approach discussed in Section 1.2. Such an exact method of solution is essential for those situations in which the true optimal solution must be determined. It is felt that many improvements could be made in the branch and bound algorithm, in addition to refinements in the computer code used to test the algorithm. It is hoped that a better method of FIXED CHARGE PROR.EM 235 determining the true maximum level of degeneracy of the convex set of solutions can be found, so that a better lower bound for the objective function at any given node may be computed. Other refinements in the handling of the degeneracy problem could probably also be found. It is suggested that research be conducted on the problem of generalizing the branch and bound approach described herein, so that it would be applicable to more general separable concave functions. Another generalization which should be considered is the case in which the fixed charges are dependent upon more than one variable. For example, a fixed charge dtj could be incurred only if the corresponding variables xt and Xj are both positive. 6. ACKNOWLEDGMENTS The author wishes to acknowledge the support of the Washington University Computing Facilities through NSF Grant G-22296 and also the Atomic Energy Commission for partial support of the research under Research Contract No. AT(11-1)-1493. REFERENCES [1] Agin, Norman, "Optimum Seeking with Branch and Bound," Management Science, 13 B- 176- 185 (Dec. 1966). [2] Balinski, M. L., "Fixed-Cost Transportation Problems," Nav. Res. Log. Quart. 8 41-54 (Jan. 1961). [3] Cooper, L. and C. Drebes, "An Approximate Solution Method for the Fixed Charge Problem," Nav. Res. Log. Quart. 14, 101-113 (Mar. 1967). [4] Cooper, L. and A. M. Olson, "Random Perturbations and the MI-MII Heuristics for the Fixed Charge Problem," Department of Applied Mathematics and Computer Science, Washington University, Report No. COO-1493-7. [5] Denzler, D. R., "An Approximative Algorithm for the Fixed Charge Problem," Nav. Res. Log. Quart. 16, 411-416 (Sept 1969). [6] Hadley, G., Linear Programming (Addison- Wesley Publishing Company, Inc., Reading, Massachusetts, 1962). [7] Hadley, G., Nonlinear and Dynamic Programming (Addison- Wesley Publishing Company, Inc., Reading, Massachusetts, 1964). [8] Robers, P. and L. Cooper, "A Study of the Fixed Charge Transportation Problem," Department of Applied Mathematics and Computer Science, Washington University, Report No. COO- 1493-9. ON COMPUTING THE EXPECTED DISCOUNTED RETURN IN A MARKOV CHAIN* David F. Hitchcock and James B. MacQueen University of California, Los Angeles ABSTRACT The discounted return associated with a finite state Markov chain X t , X 2 ■ ■ ■ is given by g(Xi) + ag(X 2 ) + a 2 g(X 3 ) + . . ., where g(x) represents the immediate return from state x. Knowing the transition matrix of the chain, it is desired to compute the expected discounted return (present worth) given the initial state. This type of problem arises in inventory theory, dynamic programming, and elsewhere. Usually the solution is approximated by solving the system of linear equations character- izing the expected return. These equations can be solved by a variety of well-known methods. This paper describes yet another method, which is a slight modification of the classical itera- tive scheme. The method gives sequences of upper and lower bounds which converge mono- tonely to the solution. Hence, the method is relatively free of error control problems. Computational experiments were conducted which suggest that for problems with a large number of states, the method is quite efficient. The amount of computation required to obtain the solution increases much slower with an increase in the number of states, N, than with the conventional methods. In fact, computational time is more nearly proportional to /V 2 , than to N :i . I. Introduction Consider a Markov chain X x , X 2 , . . ., taking values in the finite set S= [1,2, . . ., /V] that has transition probabilities p(y;x) = Pr[X t+ i — y\Xt = x]. When in state x an amount g{x) is received immediately. The discounted return associated with passing through the sequence of states X\, X2, ■ ■ -, is g{X l ) + ag(X i ) + c?g(X 3 ) + . . ., where a is the discount factor (0^a<l). The expected discounted return given Xi = x, represented by u(x) (also called the present worth of state x), is well known to be the unique solution of the equation (!) u(x) = g(x) + a^ u(y)p(y;x). An efficient machine algorithm for solving (1) is desired. The problem of computing the expected discounted return frequently arises in inventory problems, replacement problems, and elsewhere (c.f. [3], [4]). If N (the number of states in S) is reasonably small, then (1) can be solved by familiar techniques, such as Gausian elimination or the simplex method, and is a straight forward problem in machine com- putation. However, if there is a large number of states, the amount of computation, which is generally *This work was supported by the Western Management Science Institute under a grant by the Ford Foundation, and by the Office of Naval Research under Contract No. 233(75), Task No. 047-003. Computational facilities were provided by the Western Data Processing Center, University of California, Los Angeles. 237 238 D. F. HITCHCOCK AND J. B. MACQUEEN proportional to N 3 (cf.[l, p. 215]) using such methods, may be prohibitively expensive. Also, with these methods error control problems can arise even in problems of moderate size. This note describes some computational experiments using a slightly modified version of the classical iterative method for computing u. The method determines sequences of functions {u' n } and {u'n) on S which converge monotonely to u from below and above, respectively. Then u n = (u' n + u'^)/2 is taken as the estimate of u. The experiments indicate that for "large" problems (say /V ^ 100) , the amount of computation required to obtain a solution whose percentage error at each state is less than a preassigned figure, increases much more slowly than /V 3 . For example, in a number of the experiments with TV = 500, approximately 16(500 2 ) elementary (multiply and add) computations were required to achieve 0.005-percent accuracy. As with the estimate produced by the classical iteration scheme (see Eq. (5) below) convergence is found to be exponential, but at a rate significantly faster than the classical method. The monotonely convergent upper and lower bounds, which are easily calculated, provide a very satisfactory means of controlling error due to round-off, which may accumulate with other methods. The method starts with any function v\ which is at a selected state seS, with the functions t>2, t>3, . . ., being defined by (2) v„ + 1 (x) = g(x) + a 2 Vn(y)p(y; x) - (g(s) + a J v n (y)p(y; *) )• Thus v n (s) = 0, n = 1, 2, . . (3) . ... Let L' n = min x L" n = max x g(x) + a ^ v n (y)p(y; x) - v„(x) y g(x) + a ^ v n {y)p(y; x) - v„(x) The sequences {u' n } and {u" n } are then given by (4) u' n = v„+ (1- a)~ l L' n u: = v n + (l-a)-'L; The sequences {u' n } and {u'^\ have the properties: (i) u' n ss u =£ ul, (ii) u' n s£ u' n + 1 , u" n ^ u'„ +u and (hi) u' n -* u, u" n -* u. These properties are obtained directly as special cases of Theorem 2 in [2]. The modification represented by the calculation of u„ is to be compared with the classical iterative method. The latter is defined by an initial function U\ and the rule (5) u„ + i(x) = g(x) + a^ u n (y)p(y; x). If ui{x) =vi(x) + «i(s), then v n (x) = u n {x) — u n (s), and «„U) = «„(x) + V2(i-a)-'(y;+y;'), y^ = min x g(x) + a^u n (y)p(y; x) — u n {x) L v Jl=m<ix x g{x) +a^u n (y)p(y; x) —u„(x) (6) where (7) COMPUTING RETURN IN MARKOV CHAIN 239 II. Experimental Procedure To explore the computational efficiency of the suggested method, it was decided to investigate four main types of transition matrices, with TV equal to 10, 50, 100, 500, or 1,000. The discount factor was varied over the four values 0.90, 0.95, 0.99, and 0.995 representing the likely range of applications. In all the experiments, the initial function v\ was identically zero. The four main types of transition matrices were as follows: (I) Random. In this type of matrix, the row entries were randomly generated by selecting for each entry a uniformly distributed random number between and 1 and then normalizing by dividing by the row total. (II) Large probability in the left hand column. Here the first column has a specified value po (i.e., p(l; x) —po for x= 1,2, . . ., N). The value of po was varied over the three values 0.1, 0.5, and 0.9, with the remaining entries generated randomly as above (summing to 1 — po, of course). (III) Large probability on the diagonal. Here all diagonal elements were set equal to po, and the other entries were random as before. The value of po was varied as above. (IV) Large probabilities scattered. In this type of matrix, the first row was generated by taking four probabilities pi=p /2, p 2 = po/4, p 3 = po/8, and p 4 = p /8 and spacing them evenly in the first row. For example, for /V=100 we have p(l; l)=pi, p(25;l)=p 2 , p(50; l)=p 3 , and p(75;l)=p 4 . The remaining rows of p(y; x) were obtained from the first row by cyclic permutations, the last entry in the first row becoming the first entry in the second row, etc. It was decided arbitrarily to have g(x) vary in multiples of 5 as follows: for n—10, g(l) = 5, g{2) = 10, . . ., g(10) = 50. For the remaining values of N, all of which were multiples 10, the first N/10 values of g were 5, the second N/10 values were 10, etc., with the last A710 values being 50. Evi- dence that this method of choosing g, while arbitrary, was not malicious, was obtained by repeating all the experiments for N= 100 for matrices of Type III (large probabilities on the diagonal) with a different scheme for choosing g. Each value of g(x) , x— 1, 2, . . ., N, was determined as a uniformly distributed random number between and 100. Matrices of Type HI were chosen for this additional experimentation since they were the most difficult of the four in terms of convergence times. The number of iterations (see below, Table 1) using the randomly determined g was essentially identical to those obtained with the fixed g, differing by 0, 1 or at most 2 iterations. III. Results The main results are presented in Table 1. The entries in the table represent the number of iterations (the value of n in Eq. (4)) required for the maximum possible percentage error (MPPE) of the estimate u n to be less than 0.005 or 0.5 percent. Thus, for the indicated value of n, max x \u(x) — iin(x) | /u(x) is less than 0.00005 or 0.005. The upper and lower entries in each cell in the table correspond respectively to the 0.005 and 0.5 percent figures. The values of n were actually computed by allowing the computer program to run until max x {u'n{x) — u' n (x))lu' n {x) =£0.0001, which implies max x |«(x) — u(x)\lu(x) =£ 0.00005. The 0.5-percent values were then determined by inspection of the output which included the MPPE for each iteration. It should be noted that the MPPE may be considerably larger than the actual percentage error. As a typical example, consider the matrix of Type III with P = 0.90, /V= 100 and a = 0.90 where the MPPE after 16 iterations was 6.4 percent. However, ui 6 (x) differed by at most 1.4 percent from u 52 (jc), and the latter (which corresponds to the last iteration performed) has at least 0.005 percent accuracy. 240 D - F HITCHCOCK AND J. B. MACQUEEN The MPPE would not be a meaningful criterion of accuracy if u(x) were very near zero for some x. This possibility was avoided by the choice of strictly positive values of g. The MPPE can also be rendered meaningless by adding a large positive constant to each value of g, since if u is the solution with immediate return function g, u+A/(l — a) is the solution with immediate return g + A. Thus by making A sufficiently large, the MPPE in the approximate and easily computed estimate, 0. — A/(l — a) can be made arbitrarily small. Again, the choice of g avoids this trivial extreme. Several general findings were noted. First, the MPPE exhibited a strong tendency to decrease exponentially with the number of iterations, i.e., as kr n ,0 > r > 1. Moreover k is typically much less than unity. This was clearly evident on plotting the error against n on semi-logarithmic paper. The maximum percentage error for the classical method defined by equation (5), also decreases exponen- tially, approximately as Ka" where K is near unity. However, r was appreciably less than a, the ratio r/a being .9 or less in all the experiments. Thus the ratio of percentage errors for the two methods after 30 iterations, say, would be around .05 and perhaps smaller. A second finding is that the method automatically takes advantage of certain simplifying features of the transition matrices. Thus in matrices of Type I and II, the distribution of the next value of g is very similar from state to state, and the sequence g(Xi), g(X 2 ), . . ., is a sequence of similarly distributed, nearly independent random variables. For that reason, a good estimate for u(x) is simply the immediate return g( x) plus the mean value of g times a/ (I — a), and a good computational scheme ought to take advantage of this simplicity. The experimental results show that the method does. Some- what surprisingly, when the rows of the transition matrix are identical, then U2(x) =u'i(x) =g(x) +a(l-a)- 1 ^g(y)p(y; x)=u(x), y and the exact result is obtained immediately (given V\ — 0). This suggests a general rule, that the less the effect of the initial value of x on the sequence of ensuing returns, the less the amount of computation required. On this basis, one would predict the Type IV matrix to be easier to solve than the Type HI matrix, which in fact is the case. A third resul suggested by Table 1 is that the amount of computation required to achieve a given MPPE increases approximately as N 2 . There are approximately N 2 elementary computations (multiply and add being counted as an "elementary computation") per iteration, so that the amount of computa- tion required in n iterations is proportional to nN 2 . As is evident in Table 1, for a given combination of matrix type, Po and a, n increases slowly and at a decreasing rate, or does not increase at all, as /Vis increased. The worst case from the point of view of the relative increase in computation with N is the Type IV matrix with a = 0.95 and po = 0.9. However, even in this instance only nine more iterations are required to achieve 0.5-percent error as N is increased by a factor of 10 from 50 to 500. Note, how- ever, that in Table 1 with N — 10, there are a good number of instances where more than 10 iterations and hence more than N 3 elementary computations are required. Thus the relative advantage of the method as against conventional method requiring approximately /V 3 elementary computations, would appear to hold only for large values of N. Finally it is to be noted that the suggested method seems to be insensitive in an appropriate way to near zero entries in the transition matrices. This is perhaps a major reason for the rapid convergence. Type I matrices for which 0.005-percent accuracy was achieved in four iterations, particularly illustrate this result as it is obvious that for N = 500 or N = 1 ,000 the entries are extremely small. COMPUTING RETURN IN MARKOV CHAIN 241 TABLE I. Iterations Required to Achieve 0.5 and 0.005 percent Maximum Possible Error Matrix Type I II III IV \ N a \ 10 50 100 500 1000 Prob Po \ N a \ 10 50 100 500 1000 \ N \ 10 50 100 500 1000 \ N a \ 50 100 500 1000 0.995 7 4 5 4 5 3 4 3 4 3 0.9 0.995 4 3 4 3 4 3 4 3 0.995 83 46 90 51 91 51 0.995 58 30 61 32 49 32 0.99 7 4 5 4 5 3 4 3 4 3 0.99 4 3 4 3 4 3 4 3 0.99 80 45 87 49 87 49 0.99 43 19 54 24 0.95 7 4 5 4 5 3 4 3 4 3 0.95 4 3 4 3 4 3 4 3 3 3 0.95 65 37 66 38 69 39 0.95 36 16 43 20 37 25 0.90 6 4 5 4 5 3 4 3 4 3 0.90 4 3 4 3 4 3 3 3 3 3 0.90 50 28 52 30 52 29 0.90 25 12 37 17 29 20 0.5 0.995 5 3 5 3 5 3 6 3 0.995 15 9 16 9 16 9 0.995 15 9 19 10 0.99 5 3 5 3 5 3 5 3 0.99 15 9 16 9 16 9 0.99 11 6 12 7 0.95 5 3 5 3 5 3 5 3 4 3 0.95 15 8 15 9 16 9 0.95 11 6 12 7 12 7 9 0.90 5 3 5 3 5 3 4 3 4 3 0.90 15 8 15 9 14 9 0.90 10 6 11 6 12 7 8 0.1 0.995 7 4 6 3 6 3 5 3 0.995 7 4 6 4 6 4 0.995 6 4 6 3 0.99 7 4 6 3 6 3 5 3 0.99 7 4 6 4 6 4 0.99 6 4 5 3 0.95 7 4 5 3 5 3 5 3 4 3 0.95 7 4 6 4 6 4 0.95 6 4 5 3 5 3 4 0.90 7 4 5 3 5 3 4 3 3 3 0.90 7 4 6 4 6 4 0.90 5 3 5 3 5 3 4 4 4 3 REFERENCES [1] Kunz, Kaiser S., Numerical Analysis (McGraw-Hill, New York, 1957). [2] MacQueen, James B., "A Modified Dynamic Programming Method for Markovian Decision Problems," J. Math. Anal, and Appl. 14, 38-43 (1966). [3] Scarf, Herbert E., "A Survey of Analytical Techniques in Inventory Theory." In Multistage Inventory Models and Techniques, Herbert E. Scarf, Dorothy M. Gilford, and Maynard W. Shelly, editors (Stanford Univ. Press, 1963). [4] Wagner, Harvey M., Michael O'Hagen, and Bertil Lundt, "An Emperical Study of Exactly and Approximately Optimal In- ventory Policies," Management Science, 11, 690-723 (1965). MORE ADO ABOUT ECONOMIC ORDER QUANTITIES (EOQ) Victor J. Presutti, Jr. and Richard C. Trepp Operations Analysis Office Headquarters, Air Force Logistics Command ABSTRACT This paper is concerned with the determination of explicit expressions for economic order quantities and reorder levels, such that the cost of ordering and holding inventory is minimized for specific backorder constraints. Holding costs are applied either to inventory position or on-hand inventory, and the backorder constraint is considered in terms of the total number of backorders per year or the average number of backorders at any point in time. Through the substitution of a new probability density function in place of the normal p.d.f., explicit expressions are determined for the economic order quantities and the reorder points. The resulting economic order quantities are independent of all backorder constraints. It is also concluded that under certain conditions, the minimization of ordering costs and inventory holding costs (applied to inventory position), subject to a backorder constraint, is equivalent in terms of reorder levels to minimization of the safety level dollar investment subject to the same backorder constraint. INTRODUCTION We are concerned with a single-echelon, multi-item, continuous review inventory system in which the process generating demands does not change with time. When an item's reorder level is reached, an order is placed. After the elapse of a leadtime, which is assumed constant, the entire order is re- ceived. For each item, it is desired to determine the constant reorder quantity and reorder point which will minimize (summing all items in the inventory) the holding and procurement costs subject to an overall backorder constraint. The holding cost may be a function of either on-hand inventory or inven- tory position. We will constrain total units backordered or average number of units in a backorder position, which, as is well known, is also the unit-years backordered per year. In either case, the backorders may be weighted by item essentiality. MODEL DERIVATION Definitions The following definitions apply throughout the discussion: Q: order quantity (units) r: reorder level (units) D: annual demands (units) k: safety factor cr: standard deviation of unit demand during a leadtime a: holding cost factor A: ordering cost (dollars per order) c: item cost (dollars per unit) fju: mean leadtime demands Z: item essentialtiy (relative military worth) 243 244 V. J. PRESUTTI, Jr. AND R. C. TREPP x: units demanded during leadtime T)\ constraint — total essentiality-weighted units backordered /3: constraint — expected number of essentiality-weighted units in a backordered position at any point in time t: leadtime. We also require: i) r = fx + kcr ii)/(x),the probability density function (p.d.f.) describing the distribution of units demanded during the leadtime: (1) 2^ ex P ( X — /JL for all x. It is easily shown that this implies: (2) and (3) Pr[x > fx+ ko-]= 0.5 exp (- V2k) for k 2* Pr[x > fJL + k<r] = 1 - 0.5 exp ( V2A) for k < 0. The form of f(x) was inspired by the work of L. L. Parker [3] who essentially assumed demands were normally distributed, and then approximated the complementary cumulative function by a function of the form ae~ bx . We originally chose a and b such that the mean absolute deviation from the normal distribution over the range =s k =s 3 was minimized (a = 0.57, 6=1.35). We now consider, however, that a = 0.5 and b— v2 is more appropriate.* In addition to defining a p.d.f. that is symmetric about the mean (/u,), these values also cause Pr hi < ?—£ < k 2 to be independent of fi and u (a desirable attribute also shared by the normal p.d.f.). Further, with these values, the expected backorders at a random point in time for k = is the same as if demands were normally distributed, for large Q/cr. This similarity was considered desirable since we are primarily interested in average backorders (rather than total backorders) and we will often have no safety level. The value of P r \k =£ =£ k + 0.2 for k = to 4.0 is given in Table 1 for both the normal p.d.f. and f{x). We will now operate as if A.)-^-(-vi(^)r *We are indebted to Mr. Alan Kaplan of the AMC Inventory Research Office who suggested that we reconsider our values for a and b. **The absolute value of — ■ — is omitted. EOQ 245 Table 1 Comparison o/f(x) to the Normal p.d.f. K A" B" C n (xY C f (xY 0.0793 0.1232 0.0793 0.1232 0.2 0.0761 0.0928 0.1554 0.2160 0.4 0.0703 0.0700 0.2257 0.2860 0.6 0.0624 0.0527 0.2881 0.3387 0.8 0.0532 0.0397 0.3413 0.3784 1.0 0.0436 0.0299 0.3849 0.4084 1.2 0.0343 0.0226 0.4192 0.4310 1.4 0.0260 0.0170 0.4452 0.4480 1.6 0.0189 0.0128 0.4641 0.4610 1.8 0.0132 0.0097 0.4773 0.4704 2.0 0.0088 0.0073 0.4861 0.4777 2.2 0.0057 0.0055 0.4918 0.4832 2.4 0.0035 0.0041 0.4953 0.4874 2.6 0.0021 0.0031 0.4974 0.4905 2.8 0.0013 0.0024 0.4987 0.4928 3.0 0.0006 0.0018 0.4993 0.4946 3.2 0.0004 0.0013 0.4997 0.4959 3.4 0.0001 0.0010 0.4998 0.4969 3.6 0.0001 0.0008 0.4999 0.4977 3.8 0.0001 0.0006 0.5000 0.4982 4.0 0.0000 0.0004 0.5000 0.4987 — -=S/t + 0.2 \fom(x)=— i=exp •V2^ \ 2<r 2 / ■B = Pr\k* ^k+0.2\forf(x)=— c .exp(-V2 cr J 2cr \ cr x — fj. x-V- } X-/A C(x)=Pr =£ s Ar-H0.2l for n(x) and /(*), respectively. 246 V. J. PRESUTTI, Jr. AND R. C. TREPP Notice that this is strongly biased against negative k. In fact, by operating as if the probability dis- tribution could be described by Eq. (4), we are, in effect, saying that limit Pr[x > p. + kcr] = °°, a k— *— °° somewhat preposterous allegation. Nonetheless, this bias is considered desirable since in actual practice, we restrict k 3 s 0.* Expected Backorders In order to determine the expected number of units backordered per year (Bx) and the average number of units in a backorder status at a random point in time (Bt) , we need g(y) —the p.d.f. for the number of backorders at a random point in time. The following technique for obtaining g(y) was developed by Hadley and Whitin [1]. g(y)dy— probability that the number of units backordered is between y and y + dy, and £-(-Vl(«±5=£))* = the probability that the number of units backordered at time, t, is between y, y+ dy given that the inventory position (on-hand plus due-in minus due-out) is between p, p + dp at time t — t. But the probability that the inventory position is between p, p + dp is dpIQ. Therefore, ^r^e, P (-V2(^)), r * and and finally (7) g(y) = ^ exp (- V2 ( Z ^)) d " exp( -V2QI<r)). Consequently, the probability of a stockout (P ou t) is (8) | o " g(y)dy = (^|) | (1 - exp (- V2QI*)) exp (- V2k). Therefore, (9) B x = DP 0Ut = (~^j -Q- (1 - exp (- V2QI<t) ) exp (- V2k) (seeRef. [1], page 178) and (10) B T = f~ yg{y)dy = -^-^j (1-exp (- V2QI<r)) exp (- V2A). Other Inputs We must now obtain the average number of procurements per year, the expected inventory position (E[IP]), and the expected on-hand assets (E[OH]). Clearly, the average number of procure- ments is D/Q. *The nonnegativity requirement for A: is a result of the desire to make this model acceptable to practicing logisticians. EOQ 247 It is also well known that and By definition, Therefore, £[Due-in] = \x (see Ref. [1], page 187) E[IP]=r+QI2 = ix + ka + QI2. E[IP] =E[OH] +£[Due-in] -E[B/0]. E[OH] =E[IP] -£[Due-in] +E[B/0] = fi + ko- + QI2-ix + B T = k<r + QI2 + B T . FOUR INVENTORY MODELS We are now ready to proceed to the four models which will enable us to determine Q and r for each item in such a way as to minimize the holding cost (applied either to E[OH] or E[IP]) plus the pro- curement cost, subject to an overall backorder constraint. This constraint may be either the total units backordered per year (i.e., the stockout penalty is independent of the duration of the stockout) or the average number of units in a backorder position (i.e., the stockout penalty depends on the duration of the stockout). In either case, the backorders may also be weighted by item essentiality. MODEL I: Stockout penalty independent of duration of stockout; holding cost applied to on-hand inventory. We wish to minimize (LI) subject to (1.2) A AiD t - Z-nT + Z aiCi i=l v ' t=l kcri+y+^^(l-exp(-V2^) A 0.5 ZiDjO-j (l-exp(-V2^) exp (- V2ki) *S t, £V2 Q t where Z, represents item essentiality. The Method of Lagrange yields QidiCi (1.3) (1.4) *' = "V? ln -V^))(^p-XZ (J D, 2 (afim -=V2-7,,and VJ -kZiDi (1.5) o = -^ + ^i Qf 2 aidCTi V2" 1-exp [-V2 £)) 1-exp [-V2 a-iJJ V2 Qi <Ti exp (-^5) For reasons which will soon become evident, we will not solve (1.5) for Q t until later. MODEL II: Stockout penalty independent of duration of stockout; holding cost applied to in- ventory position. "AM" ( Q t Z ~q~ + Zj aiCi[fJ<i + ki(Ti + Y 248 V. J. PRESUTTI, Jr. AND R. C. TREPP We wish to minimize (III) subject to " 0.5 ZiDm I ( _ Qi\ , The Method of Lagrange yields 1 (113) "v^ ln QidjCj 0.5 (-A) Z,Z) ( fl-exp (- V2" <?.■ (H.4) A = — — — V cnatCi, and V2 Ti £i i" (115) = -^i + ^- aiCicn Vf(l-exp(-V2a)) 1 — exp -V2& MODEL III: Stockout penalty depends on duration of stockout; holding cost applied to on-hand inventory. We wish to minimize (IH.l) subject to (III.2) |)^+|;aiC i [fe(ri+^+^g(l-exp(-V2^))exp(-V2A i )] i^^(l-ex P (-V2^))ex P (-V2fe)^ i 8. The Method of Lagrange yields (III.3) kt — ;= In V2 V2(W, 0.5(7, ( 1 -exp (- V2 & ) ) (aid-kZi) (III.4) 1,(5^=^. and (III.5) °=-^ +S f- aiCicn V2(l exp V2 CTi 1 — exp V2& <?, -^expf-V^ o-j CXi MODEL IV: Stockout penalty depends on duration of stockout; holding cost applied to inventory position. EOQ 249 We wish to minimize (IV.l) subject to (IV.2) y -jt + i aiCi ( /*■• + kiai + 9 i=l v* i = 1 \ ^ The Method of Lagrange yields (IV.3) ki= ~vz ln V2 QidiCi 0.5 (- \)ZiO-, 1 -exp - V2 (J, (IV.4) (IV.5) " cr ( aiCi : > — ;^t:, and <?? 2 " atCjCTi V? (l-exp - V2" -a <Ti l-exp l-v^ 5i ^exp(- (3i cr, Notice that in all four models Qi is the same (Eqs. (1.5), (II. 5), (HI. 5), and (IV.5)). The quantity Qi is completely independent of the backorder constraint regardless of whether or not the constraint is time-weighted. We also note that we cannot solve explicitly for Q. Conceptually, there is no problem. The exact solution for Q can be obtained through the use of an iterative routine (assuming access to a computer). An alternate and perhaps a more suitable approach is to recognize that in most cases (!-- p(-Vjfi)) ] and can therefore be ignored. Consequently, AiDi (13) and (14) atCi ajCiCTi ~ ~ Qf + 2 V2C, Qt V2 + V a iCi + lv2T Notice that as <t, — > 0, Eq. (14) approaches the well known Wilson lot size formula. DISCUSSION -MODEL APPLICATION Holding Cost Considerations We wish now to discuss the difference between applying holding cost to the on-hand inventory (Models I and III) and applying it to the inventory position (Models II and IV). In spite of tradition, we feel it is more logical to apply the holding cost to the inventory position (see Ref. [4] for an argument supporting this point of view). From a practical point of view, the difference in results is small; however, the inventory position models have two very desirable attributes. First, there is an explicit solution *It is easily shown that this value is an upper bound for the exact solution. 250 V. J. PRESUTTI, Jr. AND R. C. TREPP for X, the Lagrangian, in these models — which is not the case when the holding cost is applied to the on-hand inventory. Second, in the inventory position models the safety level expressions are equivalent to those that would have been obtained had we somehow arrived at the buy quantities ((Vs) and then independently minimized the cost of holding safety stock subject to the same backorder constraints used in Models I through IV. That is, under the inventory position models, the inventory manager who minimizes backorders subject to a safety level cost constraint would arrive at A:,'s equal to the A;,'s that he would have calculated had he minimized the total costs of ordering and holding inventory subject to a backorder constraint. This is especially important when the optimal (Vs are not used: for example, when the policy is to purchase 3, 6, 9, or 12 months' worth of stock, or even when Qi is determined by the Wilson lot size formula. If the same holding cost factor is applied to all items (i.e., at = a for all t), the results obtained in models II and IV are equivalent to those determined when minimizing the safety level investment subject to a backorder constraint. Negative k's Although all four models are strongly biased against negative &'s, as long as there is a heterogeneous inventory with some high cost, low demand items, there will be some items for which the models will specify negative &'s. It could be argued that if even with a strong bias negative A:'s are obtained they should be accepted. We have no philosophical quarrel with this argument; however, in this paper, we insist that there be no negative safety levels. In order to accomplish this, we will set all negative A's equal to zero. Consequently, we will have more protection than the models call for and thus have fewer backorders, but a higher cost. The problem is to determine the X that will allow, taking into account the nonnegativity of k, an acceptable number of backorders at a reasonable cost. Since we believe that Model IV is most appropriate for the Air Force's nonreparable items, we will use it to describe the methodology for determining A. Employing a simulation using actual data for about 5,000 Air Force items (for which seven years of history is available) we will select a range of values for [—A.]. Each value of [—A] will yield 5,000 A;'s — one for each item. The nonnegative A:'s will be accepted and the negative &'s will be set equal to zero. We will then, using these A;'s which are now all nonnegative, obtain the average number of essentiality-weighted units in a backorder status divided by the number of essentiality-weighted units demanded annually (call this ratio /?«), and the cost of safety stock divided by the dollar value of annual demands (call this ratio R t ). Subsequently, Figures 1 and 2 will be plotted as shown below: -3t Given Figures 1 and 2, the inventory managers can select either R B or R$ (hopefully they will make the selection by considering both the cost and the effectiveness aspects) and then choose from the appropriate graph the value of [— X]. The ability to choose on the basis of fund availability is extremely important since, in practice, the inventory manager is considerably more likely to be subject to a system constraint of dollars available for safety level investment (out of pocket cost) than to a back- eoq 251 order constraint. Hadley and Whitin emphasize this point in their* discussion of budget constraints. They conclude: "None of the inventory models discussed in operations research literature appears to provide an adequate treatment of budgetary constraint. In the first place, budgets are usually set in terms of actual dollars to be spent, whereas the costs included in inventory control models are often not the out-of-pocket cost type."* An example of how the inventory manager (IM) might use these graphs (possibly to reduce "actual dollars to be spent") is as follows: Suppose the IM submits a budget request, for operating and safety stock, which allows 20 percent of the dollar value of annual demands for safety stock. If he receives the monies requested he can obtain from Figure 2, the [— A.] that will yield an Rt of 20 percent. If he gets M dollars less than requested, and that M dollars represents 5 percent of annual dollar demands, he can set R% equal to 15 percent (20% — 5%) and obtain the new value for [— \]. While this paper doesn't concern itself (other than the portion immediately preceding) with constrained optimization, we have applied the demand distribution described in this paper to models that do constrain the buy quantities (i.e., the (?'s), the on-hand inventory, and the inventory position as well as the safety stock. None of these methods, however, is satisfactory when it comes to reducing actual dollars to be spent by very large amounts below computed requirements. All of these constraining models reduce dollars spent by reducing the reorder level (i.e.. by reducing k) and/or by reducing the order quantities in some manner. We cannot, however, reduce the reorder level greatly without paying a very large price in terms of backorders. If most of our expenditures were for items with multi-year buy quantities, then we could temporarily reduce actual dollars to be spent by reducing these buys to a year or less. Unfortunately, by far the greatest majority of dollars are spent on items for which the buy quantity is on the order of three months worth of stock — which theoretically eliminates the pos- sibility of obtaining a large reduction in actual dollars to be spent by reducing the buy quantity. Backorder Costs Known The reader will note that in the event the stockout cost (s,) is known, the problem is simplified. The (?,'s remain unchanged and the A,'s of each of the models have a 5; in place of [— A] Z, (see Eqs. (1.3), (II.3), (III.3), and (IV.3)). There is no [- A] . Final Note The Air Force is planning to test all four models, nevertheless, at the present time, it favors (un- officially) Model IV. This model has also been brought to the attention of the Department of Defense by the Air Force for consideration during a series of recent joint services meetings designed to further standardize inventory policies and objectives among all DoD components. REFERENCES [1] Hadley, G. and T. M. Whitin, Analysis of Inventory Systems (Prentice Hall, Englewood Cliffs, New Jersey, 1963). [2] Hadley, G. and T. M. Whitin, "A Review of Alternative Approaches to Inventory Theory," The RAND Corporation, Santa Monica. California. RM-4185-PR (Sept. 1964). (3| Parker, L. L.. "Economical Reorder Quantities and Reorder Points with Uncertain Demand," Nav. Res. Log. Quart.. //, 351-358 (Dec. 1964). [4] Schrady, D. A., "Continuous Review Inventory Practice in the Navy Application," Report #NPS55o9031A, United States Naval Postgraduate School, Monterey, California (Mar. 1969). * Reference [2|, page 152. THE RELATIONSHIP BETWEEN DECISION VARIABLES AND PENALTY COST PARAMETER IN (Q, R) INVENTORY MODELS Alan J. Kaplan AMC Inventory Research Office Army Logistics Management Center Fort Lee, Virginia ABSTRACT This paper is concerned with the optimum decision variables found using order quantity, reorder point (Q, R) inventory models. It examines whether the optimum variables (Q* and R*) are necessarily monotonic functions of the backorder cost parameter (or equivalently of the performance objective). For a general class of models it is proved that R* must increase as the performance objective is raised, and an inequality condition is derived which governs how Q* will change. Probability distributions of lead time demand are cited or found for which Q* increases, Q* decreases, and Q* is independent of increases in performance objectives or backorder cost parameter. INTRODUCTION This paper is concerned with the optimum decision variables found by using order quantity, reorder point (Q, R) inventory models. It examines whether they are necessarily monotonic functions of the penalty (backorder) cost parameter (or equivalently of the performance objective). It attempts to make general statements about the direction of such monotonicity. Such information is of particular interest in contexts where the true cost of a backorder is difficult to estimate — and the performance objective is determined in large part by its impact on inventory levels. It is easy to be misled by intuition. Our own empirical experience was with a (Q,R) model in which backorder cost is proportional to the length of time a backorder is outstanding, and lead time demand is normally distributed. For this case, the optimum reorder point (R*) always increased as the penalty cost parameter (P) increased, while the optimum order quantity (Q*) decreased. This is illustrated in Herron's [2] dimensionless graphs for the normal distribution. The conclusion can be proved analyti- cally for deterministic demand and penalty proportional to time on backorder by solving directly for Q *; e.g., see Hadley and Whitin's [1] solution. On the other hand, Victor Presutti [3], found Q* to be independent of the penalty cost parameter when he used a pseudo-exponential distribution for lead time demand. t That finding prompted this paper. We initially assume penalty cost is proportional to the length of time a backorder is outstanding and begin by reviewing the relevant aspects of the appropriate (Q, R) model. Next we obtain some useful theoretical conclusions, including a proof that R * must always increase when P does and an inequality condition which governs how Q* will behave. We apply these results, and then in the last section extend them to other {Q, R) models. t/(x) = A.ie x " forx^O, and/(jc) denned elsewhere so that / f(x)dx=\. Presutti's work is based on Lawrence Parker's approximation to the normal complementary cumulative. 253 254 A - J- KAPLAN The theoretical conclusions require the use of an approximation for average backorders (Eq. (2)). This approximation is commonly used in the actual computation of /?* and Q* [1], [2J, and does not negate the value of the results. THE MODEL The model is the basic steady state, continuous review (Q, R) model applicable when penalty cost is proportional to the length of time a backorder. is outstanding. For the model to give strictly optimal answers in the general case, procurement lead times must be constant and demands must be for one unit at a time. However, the model is being used profitably even when these conditions do not hold. For ease of exposition we will assume that the lead time demand distribution is continuous. One good reference, particularly for providing the basis for Eqs. (1) and (2), is Hadley and Whitin [1]. Let TC = average total cost per year; H= holding cost per year per unit of stock on hand; P= penalty cost per year per unit of stock on backorder; A = administrative and/or setup cost; A — average yearly demand; A/. = expected demand during procurement lead time; B — average amount on backorder at a random point in time; G(x) — probability that lead time demand equals or exceeds x; L = lead time. Then TC(R,Q)=A(klQ)+H- (R + QI2-\,. + B) +P • B = A • {KIQ)+H- (R + QI2-k,.)+ P'B (1) P' = P + H B = ±r (x-R)[G(x)-G(x + Q)]dx. v Jk For computational purposes, the following approximation is often used for B: (2 ) B = jjf~(x-R)G(x)dx. The approximation is valid to the extent that the probability of demand exceeding R + Q is small. We will utilize this expression for B. THEORETICAL RESULTS Let Qi, R, be the optimal decision variables for P'=/ > t ; Qi, Ri be the optimal decision variables for P' =P 2 ; Pz>Pr, Bt be expected average backorders for decision variables Qi, Rr, F(Q,R) (Q,R) INVENTORY MODELS 255 dB/dQ SB/dR' The following two theorems will be proved: THEOREM 1: R* (the optimal R) must increase as P increases. THEOREM 2: The sign ofdF(Q, R)/dR, if it is consistent, determines whether ()* will also increase. The key to these results is the corollary to Lemma 1 that the algebraic sign of Q 2 — Q\ is determined by the algebraic sign of F(Q 2 , R 2 ) — F(Q U R t ). Lemma 2 and its corollary demonstrate the intuitively "obvious": that expected backorders (B) must decrease in the optimum as P increases and, therefore, that it is not possible for both Q* and R* to decrease. Lemma 3 shows that dF{Q, R)ldQ is always negative. Lemma 2 is used in proving Theorem 1, while Lemma 3 is used in proving Theorem 2. LEMMA 1: The algebraic sign of F(Q 2 , Rz) — F{Q\, R\) is the same as the algebraic sign of Qi~Q\- (A zero value for one implies a zero value for the other.) PROOF: Taking the partial derivatives of the total cost equation, Eq. (1), with respect to Q* and /?*, and setting the derivatives to zero gives two conditions: (3a) (dBldQ*)-P'-^-AI(Q*y + H/2 = 0; (3b) (dB/dR*)-P' + H=0. By algebraic manipulation (transposing and dividing one equation by the other), this gives F(Q*,R*) = [A -AI(Q*) 2 -HI2] +-H. Hence F(Q 2 , R 2 ) -F(Q U /?,) = [X-AKQ.r-X-AKQ.y] --//. If Q 2 exceeds Qi(i.e., Q2 — Q1 is positive), then both the numerator and denominator on the right hand side of Eq. (4) are negative; F(Q 2 , R 2 ) — F(Qi, Ri) is then positive. By similar reasoning the rest of Lemma 1 is shown to hold. COROLLARY: The algebraic sign of Qt — Q\ is the same as the algebraic sign of F(Q 2 , R 2 ) — F(Q\, Ri). A zero value for one implies a zero value for the other. LEMMA 2: The expected value of backorders (in the optimum) must decrease as P increases. PROOF: To simplify the explanation let Oi stand for operating costs corresponding to decision variables (Q u Ri) so that Total Cost = O t + B t ■ P' . We use proof by contradiction. We initially assume B 2 > B t . By the optimality of Q u Ri for/ > '=P 1 , Using P, > P, and B>>Bi, we find O. + BrP^BriPi-Pt) > O. + fi.P. + fi.^Pi-P,). Hence 2 + B 2 P 2 >0 1 + Bi P 2 . But this last inequality contradicts the optimality of Q 2 , R 2 when P' =P 2 . COROLLARY: Both Q * and R * cannot decrease as P increases (proved by Lemma 2 and inspection of Eq. (2)). LEMMA 3: dF{Q, R)/dQ is always negative. PROOF: Using Eq. (2), we get by calculus, 256 A. J. KAPLAN (5a) % = ~7p\l ('-WW** dB 1 f- v , = ~Q in G(x)dX - and (5b) Hence, dR (6) F(Q, R) ~j^(x-R)G(x)dx ±j^ (x-R)G(x)dx g — I™ G{x)dx !" G(x)dx [ X G{x)dx (JjR JR JR The denominator does not depend on Q. Clearly, the numerator decreases as Q increases. THEOREM 1: R* increases as P increases. PROOF: Using Eq. 6, we have (7a) F(Q 2 ,R 2 )=B 2 +\G(x)dx J «2 and (7b) F(Q u R 1 )=Bi^r('G(x)dx. JRi Suppose R 2 < Ri. Then the denominator of (7a) will be larger than the denominator of (7b). By Lemma 2, B 2 <B\, so this means F{Q 2 , R 2 ) <F(Qi, Ri). By the corollary to Lemma 1 this implies Qz<Q\. By the corollary to Lemma 2, however, it is not possible for Q 2 < Q\ if R 2 < Ri since P 2 > Pi. Proof by contradiction. THEOREM 2A: If dF(Q, R)/dR is consistently negative, then Q* must decrease as P increases. Since R 2 >Ru and dF(Q, R)/dR is given as negative, F(Q U R 2 ) <F(Q l , R^. Now if Q 2 ><?,, then by Lemma 3 this means F(Q 2 ,R 2 ) < F(Qi, R 2 ). By combining the two inequalities, we get F(Q 2 , R 2 ) < F(Qi, Ri). However, by the Corollary to Lemma 1, this implies Q\>Q 2 . Proof by contradiction. THEOREM 2B: If dF(Q, R)/dR is consistently positive, then Q* increases with P. The proof is analagous to that of Theorem 2 A. Since R 2 > Ri, and dF(Q, R)ldR is given as positive, F(Q U R 2 ) >F(Q U /?,). Now if Q 2 <Q U then by Lemma 3 this means F{Q 2 ,R 2 ) >F(Q U R 2 ). By combining the two inequalities we get F(Q 2 , R 2 ) > F(Q t , R\). But by the Corollary to Lemma 1 this implies @i < (?2. Proof by contradiction. THEOREM 2C: If F{Q, R) is independent of R, then Q* is independent of P. For this case F(Q U R 2 )=F(Qi, Rt). If <? 2 < <?i then by Lemma 3, F((? 2 , /? 2 ) > F{Q X ,/?,), which by the Corollary to Lemma 1 implies Q 2 >Q\. Similar reasoning shows that the assumption ()i < Q 2 also leads to a contradiction. Applications There is a strong temptation to reason intuitively that Q* must increase as penalty cost increases: higher reorder quantities reduce the probability of stock out; this advantage becomes more significant as the penalty cost parameter is increased, while the impact of Q on ordering and holding costs is unchanged. Why is this reasoning wrong? The flaw is that when penalty cost increases, it becomes worthwhile to increase /?*, so that chang- ing Q* has less impact on backorders than before, and therefore may have less impact on penalty (Q,R) INVENTORY MODELS 257 costs (the product of backorders and the penalty cost parameter). The decisive factor relates to the relative effectiveness in reducing backorders of augmenting Q* vs augmenting R*. As R increases, does the relative effectiveness of augmenting Q increase? Only then will Q* increase as P increases. It is difficult to apply Theorem 2 analytically to most probability distributions, i.e., to determine analytically the sign (+, — , or mixed) of dF(Q, R)dR. No attempt was made to carry on an extensive computer aided mapping of sign vs probability distribution We did attempt to find a case when Q* would increase as penalty cost increases. Recall that for the normal distribution Q* decreases and that for exponential type distributions Q* is unchanged. Let the Weibull distribution be defined: f(x) = aXx"- 1 exp (- Xx a ); x^0,Xa^0. Then for a= V2, we have analytically shown that dF(Q, R)/dR is positive for R 3= 0, i.e., this is a case where Q* increases with R. Empirical work indicates that this conclusion about the Weibull distribution also holds true when the precise expression for average backorders (Eq. (1)) is used.t For a=l the Weibull reduces to the exponential. It is simple to show dF(Q, R)/dR — for the exponential regardless of which expression for B is used. As a matter of interest, the Weibull with a= V2 was originally examined because of an inference made from an analogy between determination of the sign of dF(Q, R)/dR, and the sign of the deriva- tive of the failure rate function of reliability theory, the latter being a standard problem. EXTENSIONS The general linear form for penalty cost is P + P • t where t is the time on backorder. This trans- lates into total expected penalty cost =P ■ B + P ■ B or, in words (cost per unit backordered per year on backorder) X (average amount on backorder) + (fixed cost per backorder occasion) X (number of demands backordered each year). That part of the total cost equation, Eq. (1), dependent on B or B becomes P' B + P ■ B where P' = P + H. The theory derived in this paper still holds provided that as the penalty cost parameters are increased, P remains equal to k • P' for constant k (or P — 0, the special case already considered). If the cost item H • B is left out of the total cost equation, Eq. (1), the condition is simply that P — k • P, which is more natural. H • B is normally very small, and in fact does not appear as a cost if holding cost depends on inventory position (including net due-in) rather than on hand. An approximate expression must be used for B as well as B. 1 1 Let fjl'W B T =kB + B. For F T {Q, R), defined in terms of dBrldQ and BB T lBR, we have Br Ft(Q,R) = (*G(x)dx + k-k-G(R) tMiss Chung-Mei Ho did much of the analytical investigation on the application of Theorem 2, and Mr. Robert Deemer did the empirical work for the Weibull. ttHadley and Whitin, chapter 4 has a derivation. 258 A. J. KAPLAN By inspection, all three lemmas proved hold with the substitutions of Bt for B, P' • Bt for P'B, and F T (Q, R) for F(Q, R). Therefore, the two theorems hold with these substitutions. Using the precise expression for B causes difficulties even when P is 0. Then B-T [y-(R + Q)]g(y)dy , C( . F{Q,R)=—±** , where g ( y ) =-&&-. [G(x)-G(x + Q)]dx . y It is no longer apparent that Lemma 3 must always hold, or that in the proof of Theorem 1 that the numeraton F(Q 2 , Rz) < numerator a£F{Qu Ri) if Qi > Q\- REFERENCES [1] Hadley, G. and T. M. Whitin, Analysis of Inventory Systems (Prentice Hall, Englewood Cliffs, N.J., 1963). [2] Herron, David P., "Use of Dimensionless Ratios to Determine Minimum-Cost Inventory Quantities," Nav. Res. Log. Quart. 13, 167-176(1966). [3] Presutti, V. J. and R. C. Trepp "More Ado about Economic Order Quantities (EOQ)", Nav. Res. Log. Quart. 17, 242-250 (1970). NEWS AND MEMORANDA PROFESSIONAL DESIGNATION IN LOGISTICS MANAGEMENT The School of Systems and Logistics, Air Force Institute of Technology, in cooperation with the Society of Logistics Engineers (SOLE), has announced the beginning of a program of professional designation for qualified military and Civil Service logistics managers of the Department of Defense. Those who meet the criteria may be awarded an appropriate certificate of "Professional Designation in Logistics Management." Certification is not restricted only to members of the Society. The purpose of the program is to recognize the efforts of the individual to increase his educa- tional background in the broad world of logistics. It recognizes that true logistics management capa- bility must encompass more than expertise in one of the several disciplines which comprise the logistics spectrum. The idea is that logistics is inherently an integrating process of these several disciplines and that the effective logistics manager will be one who has knowledge in, and understands the view- point of, the various sub-functions which the separate disciplines represent. The professional stature of the functional specialists is recognized, but this program is aimed at the need lor a portion of these managers to become logistics generalists to perform the essential function of integrating the efforts of many functional specialists. Qualification for certification rests upon the types and forms of education and experience which contribute to the overall concept of integrative management in logistics. The courses of study outlined in the qualification listing were tailored with this concept in mind. The objective is to assure that those who qualify have progressed from one of the functional specialties to a broader realization of logistics management needs through exposure to education in related and supporting disciplines and functional areas. Thus, those who qualify for certification will have demonstrated a high degree of initiative in acquiring a broad educational and experience background. They will be considered to have demon- strated their competency for advancement to higher levels of responsibility in logistics management where their broadened point of view would aid decision-making of more cost-effective nature. The certification program provides for two levels of certificate award. The first is "Professional Designation in Logistics Management" and the second is "Professional Designation in Logistics Management — Advanced." The qualification criteria includes courses in five categories as indicated: \ Certificate Advanced Category A Logistics Integration 2 2 Category Bl Materiel Acquisition 1 2 Category B2 Materiel Distribution 2 2 Category B3 Materiel Maintenance 1 2 Category C Management Techniques 2 2 Total Required 8 10 Some allowance is made in the program for waiving certain educational course requirements on the basis of field experience. Such waiver, for a maximum of two courses, must be substantiated by the applicant. The determination of whether or not to accept the experience equivalency is left to the discretion of the certification committee of the School of Systems and Logistics. Applicants are urged 259 260 NEWS AND MEMORANDA to explain their request for equivalency and waiver by an adequate explanation of the experience (where, for how long, in what job, at what level, and so on). It is also likely that many will have taken higher level educational courses than those listed in the certificate categories. These may be accepted in lieu of the described courses if the applicant includes proof of completion along with a brief outline of material covered. Certificates will be approved by the School of Systems and Logistics and the SOLE National Education Committee through a certification board process. Applicants, who need not be members of SOLE but who must be members of DOD, are encouraged to submit their requests for certification to the School of Systems and Logistics (AFIT-SL), W-PAFB, Ohio 45433. Each request must be accompanied by the necessary back-up data. CHANGE IN OFFICE OF NAVAL RESEARCH ADDRESS Effective with the relocation of Department of the Navy activities from the Main Navy Building in Washington, the address for the Office of Naval Research, including the Naval Research Logistics Quarterly, is Arlington, Virginia. Correspondence for the Naval Research Logistics Quarterly should be addressed as follows: Managing Editor Naval Research Logistics Quarterly Office of Naval Research Arlington, Virginia 22217 ■fr U.S. GOVERNMENT PRINTING OFFICE i 1970 OL-386-049 INFORMATION FOR CONTRIBUTORS The NAVAL RESEARCH LOGISTICS QUARTERLY is devoted to the dissemination of scientific information in logistics and will publish research and expository papers, including those in certain areas of mathematics, statistics, and economics, relevant to the over-all effort to improve the efficiency and effectiveness of logistics operations. Manuscripts and other items for publication should be sent to The Managing Editor, NAVAL RESEARCH LOGISTICS QUARTERLY, Office of Naval Research, Arlington, Va. 22217. Each manuscript which is considered to be suitable material tor the QUARTERLY is sent to one or more referees. Manuscripts submitted for publication should be typewritten, double-spaced, and the author should retain a copy. Refereeing may be expedited if an extra copy of the manuscript is submitted with the original. A short abstract (not over 400 words) should accompany each manuscript. This will appear at the head of the published paper in the QUARTERLY. There is no authorization for compensation to authors for papers which have been accepted for publication. Authors will receive 250 reprints of their published papers. Readers are invited to submit to the Managing Editor items of general interest in the field of logistics, for possible publication in the NEWS AND MEMORANDA or NOTES sections of the QUARTERLY. NAVAL RESEARCH JUNE 1970 LOGISTICS VOL. 17, NO. 2 QUARTERLY NAVSO P-1278 CONTENTS ARTICLES Page Markov Chain Analysis of a Situation Where Cannibalization is the Only Repair Activity by A. J. Rolre Some Estimates of Reliability Using Interference Theory by M. Mazumdar The Optimal Burn-In Testing of Repairable Equipment by J. M. Cozzolino Goal Programming in Econometrics by W. A. Spivey and H. Tamura i Production and Employment Scheduling in Multistage Production Systems by C. Haehling von Lanzenauer A Mutual Primal-Dual Linear Programming Algorithm by M. Y. Harris Interception in a Network by R. D. Wollmer The Fixed Charge Problem by D. I. Steinberg On Computing the Expected Discounted Return in a Markov Chain by D. F. Hitchcock and J. B. MacQueen More Ado About Economic Order Quantities (EOQ) by V. J. Presutti, Jr. and R. C. Trepp The Relationship Between Decision Variables and Penalty Cost Parameter in (Q,R) Inventory Models by A J. Kaplan News and Memoranda OFFICE OF NAVAL RESEARCH Arlington, Va. 22217